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0.1  Objectives 

The  objective  of  this  research  was  to  develop  approach  to  control  of  mixing  in  shear  flows  and  apply  it 
to  control  of  pattern  factor  and  thermoacoustic  instabilities  in  military  aeroengines.  A  more  general 
objective  was  to  develop  tools  for  modeling,  analysis,  and  control  design  for  unsteady  non-equilibrium 
flow  phenomena  relevant  to  operation  of  aeroengines. 

0.2  Summary  of  Accomplishments 

The  original  objective  of  developing  methods  for  control  of  mixing  relevant  to  pattern  factor  control  was 
accomplished.  Methods  for  enhancing  mixing  in  jets  in  cross  flow  using  flow  control  were  developed 
and  demonstrated  in  experiments.  Control  of  Pattern  Factor  using  jet  in  cross  flow  control  was 
demonstrated  in  a  rig  experiment.  This  work  is  described  in  Section  2.1. 

The  objective  of  controlling  thermoacoustic  instabilities  via  mixing  control  was  not  accomplished. 
The  funding  was  re-directed  towards  more  promising  fuel  control  of  thermoacoustic  instabilities.  The 
decision  was  based  on  the  fact  that  the  fuel  flow  control  was  considered  more  feasible  than  air  flow 
control  and  more  likely  to  be  transitioned  to  a  6.2  project  on  Active  Screech  Control  (jointly  funded 
by  Air  Force  Research  Lab  and  Pratt  &  Whitney). 

In  the  area  of  fuel  control  severed  accomplishements  are  worth  mentioning.  First,  a  hierarchy  of 
models  for  control  of  thermoacoustic  instabilities  was  developed.  A  method  of  control  of  flame  front 
was  demonstrated  in  a  distributed  model.  Control  of  rotating  waves  arising  as  results  of  thermoa¬ 
coustic  instability  on  a  annular  domain  was  demonstrated  in  a  reduced  order  model.  The  effect  of 
combustion  on  the  fluid  dynamics  was  analyzed  in  a  distributed  and  reduced  order  models.  The  results 
are  presented  in  Section  2.2.  Other  accomplishments  include  analysis  of  impact  of  symmetry-breaking 
(Section  2.3)  and  external  noise  (Section  2.4)  on  thermoacoustic  instabilities.  The  study  of  the  funda¬ 
mental  limitations  of  achievable  performance  described  in  the  next  paragraph  was  motivated  by  the 
control  of  thermoacoustics. 

The  most  important  theoretical  accomplishment  of  the  current  research  was  establishment  of  a 
framework  for  studying  fundamental  limitations  of  achievable  performance  in  control  of  oscillations 
described  by  nonlinear  models,  including  delays,  and  driven  by  broad-band  noise  described  in  Sec¬ 
tion  1.1.  The  framework  is  based  on  frequency  domain  formulation  of  model  response.  While  linear 
dynamic  components  (oscillators  and  delays)  are  easy  to  handle  in  the  frequency  domain,  the  challenge 
was  the  treatment  of  static  nonlinearities.  This  was  accomplished  by  replacing  the  nonlinearities  by 
their  Random  Input  Describing  Functions.  This  approach  was  very  effective  in  studying  limitations 
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of  performance  in  fuel  control  of  thermoacoustic  instabilities  using  on-off  fuel  valves.  The  framework 
involved  approximations.  As  a  first  step  towards  a  fully  rigorous  framework  for  analysis  of  nonlinear 
ocillations,  the  Spectral  Balance  approach  (Section  1.2)  was  introduced  and  demonstrated  in  an  ex¬ 
ample  of  a  nonlinear  model  with  multiple  attractors.  Other  accomplishments  include  introduction  of 
two  methods  of  adaptive  control  of  oscillations  with  uncertain  parameters  (Section  1.3)  and  analysis 
of  uncertainty  propagation  in  complex,  interconnected  dynamical  systems  (Section  1.4). 

Finally,  a  linear  framework  for  control  of  wave  phenomena  on  annular  domain  was  established. 
The  flutter  control  problem  described  in  Section  2.5  was  used  to  motivate  the  study. 

The  results  of  the  current  research  are  summarized  in  11  journal  papers  (5  published,  3  in  print,  3 
submitted)  and  20  conference  papers.  Three  invited  sessions  were  organized  with  support  from  current 
grant. 


0.3  Organization  of  the  report 

We  provide  extended  abstract  for  the  results  that  are  published  and  hence  easily  available.  We  provide 
more  details  for  the  results  that  are  not  available,  like  papers  that  are  submitted  or  in  the  case  of 
ungoing  research.  List  of  references  is  available  in  the  last  two  chapters  of  the  report.  References  to 
journal  papers  that  were  written  under  this  contract  are  referenced  with  letter  “j”  in  front  ([jl],  [j2], 
etc.).  Conference  papers  written  under  this  contract  are  referenced  with  letter  “c”  in  front  ([cl],  [c2], 
etc.). 
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Chapter  1 

Control  Theory  for  Nonlinear 
Oscillations  in  Military  Aeroengines 

1.1  Analysis  of  fundmental  limitations  of  achievable  performance  in 
control  of  oscillations  in  nonlinear  systems 

The  most  important  theoretical  accomplishment  of  the  current  research  was  establishment  of  a  frame¬ 
work  for  studying  fundamental  limitations  of  achievable  performance  in  control  of  oscillations  described 
by  nonlinear  models,  including  delays,  and  driven  by  broad-band  noise.  The  framework  is  based  on 
frequency  domain  formulation  of  model  response.  While  linear  dynamic  components  (oscillators  and 
delays)  are  easy  to  handle  in  the  frequency  domain,  the  challenge  is  the  treatment  of  static  nonlin¬ 
earities.  This  was  accomplished  by  replacing  the  nonlinearities  by  their  Random  Input  Describing 
Functions .  This  approach  was  very  effective  in  describing  limitations  of  performance  in  fuel  control 
of  thermoacoustic  instabilities  using  on-off  fuel  valves.  The  framework  involved  approximations.  As 
a  first  step  towards  a  fully  rigorous  framework  for  analysis  of  nonlinear  ocillations,  the  Spectral  Bal¬ 
ance  approach  was  introduced  and  demonstrated  on  an  example  of  a  nonlinear  model  with  multiple 
attractors. 

1.1.1  Fundamental  limitations  of  achievable  performance  in  control  of  thermoa¬ 
coustic  instabilities 

Thermoacoustic  instabilities  in  gas  turbine  and  rocket  engines  develop  when  acoustic  waves  in  com¬ 
bustors  couple  with  an  unsteady  heat  release  field  in  a  positive  feedback  loop.  Fuel  control  was 
demonstrated  to  be  an  effective  way  of  reducing  the  level  of  pressure  oscillation  in  combustors  [46]  [6] 
[7]  [30]  [23]  [17]  [11]  [24]  [45].  However,  the  achieved  reduction  of  pressure  oscillation  between  experi¬ 
ments  ranges  from  6dB  to  20dB  [11]  [24].  Moreover,  in  some  cases  the  attenutation  of  the  oscillation 
at  primary  frequency  is  accompanied  by  excitation  of  the  oscillation  in  some  other  frequency  band  [6] 
[7]  [30]  [17]  [45].  This  phenomenon  is  refered  to  as  secondary  peaks  or  peak  splitting. 

Until  recently,  the  question  of  what  are  the  factors  that  impact  the  achievable  level  of  attenuation 
of  pressure  oscillations  with  active  fuel  control  was  not  addressed.  In  particular,  the  cause  of  the  peak 
splitting  phenomenon  (two  peaks  in  the  pressure  spectrum  with  control)  observed  in  many  combustion 
experiments  [6]  [7]  [30]  [17]  [45]  did  not  have  a  satisfactory  explanation.  One  of  the  reasons  was  that 
the  majority  in  the  thermoacoustic  instability  control  community  believed  that  the  thermoacoustic 
instability  arises  only  as  a  stability  loss  leading  to  a  limit  cycling  behavior  (one  peak  in  the  pressure 
spectrum),  and  any  control  action  stabilizing  the  equilibrium  of  the  thermoacoustic  system  will  result 
in  quenching  the  oscillations  (no  peaks). 
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Experimental  results 


Simulation 


Figure  1.1:  PSD  of  pressure  signal  with  on/off  control  of  one,  two,  or  three  liquid  fuel  nozzles.  The 
peak  splitting  phenomenon  observed  in  experiment  reproduced  in  reduced  order  model  simulation. 


Two  major  contributions  of  the  UTRC  combustion  control  team  (including  academic  partners) 
were: 

1.  A  method  for  determining  whether  a  given  combustor  should  be  modeled  as  a  limit  cycling  or  a 
noise-driven  stable  system  [37,  38]. 

2.  Analysis  of  factors  that  affect  the  achievable  level  of  suppression  of  pressure  oscillations  for  both 
stable  and  limit-cycling  systems  [13,  12,  3,  4,  2]. 

The  analysis  of  the  data  from  combustion  experiments  has  shown  that  industrial  combustors  are  often 
described  by  stable,  linear,  noise-driven  models.  Standard  frequency  domain  methods  were  applied 
to  such  models  showing  that  in  some  conditions  the  transport  delays  and  limited  actuator  bandwidth 
fundamentally  limit  the  level  of  achievable  suppression  of  the  pressure  oscillations  [3].  The  origin  of 
the  peak-splitting  was  explained  using  stable  linear  models  with  large  delay  in  actuation  due  to  the 
fuel  transport  process. 

Under  the  current  contract,  the  frequency  domain  fundamental  limitation  were  extended  to  non¬ 
linear  combustion  models,  including  the  limit  cycling  ones  [4,  2]. 

The  peak-splitting  phenomenon  was  observed  in  DARPA-sponsored  sector  rig  experiments  and  was 
reconstructed  in  a  model  (see  Figure  1.1).  Because  on/off  valves  were  used  for  control,  a  linear  analysis 
was  not  applicable.  However,  a  nonlinear  analysis  involving  Random  Input  Describing  Functions 
[13,  12,  2,  24,  3,  4]  allowed  an  explanation  of  the  phenomenon.  It  has  been  shown  that,  as  in  the  linear 
model  case,  the  pressure  oscillations  cannot  be  arbitrarily  suppressed  due  to  non-minimum  phase 
effects  (transport  delay)  and  limited  actuator  bandwidth.  It  was  also  shown  that  the  effects  of  the 
driving  disturbances  and  saturation  nonlinearities  needs  to  be  incorporated  in  the  analysis.  The  sector 
rig  model  was  derived  in  the  form  of  a  feedback  interconnection  of  a  stable  linear  transfer  function 
Go(juj),  and  a  relay  nonlinearity  /(•),  subject  to  a  driving  disturbance  N(jcj)  as 

X(jw)  =  G0(jcu)(N(juj)  -  Y(ju))  (1.1) 

y(t)  =  f(x{t)).  (1.2) 

for  Go {juj)  =  Gi(ju))Gc{joj)-  Gi(jcu)  represented  the  serial  connection  of  the  thermoacoustic  response 
to  fuel  valve  output  dominated  by  a  single  resonant  mode  with  a  natural  frequency  of  about  200Hz 
and  a  non-minimum  phase  transfer  function  representing  the  transport  delay  from  the  fuel  injection 
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location  to  the  flame.  Gc{ju>)  wa s  the  transfer  function  of  a  phase-shifting  controller  used  in  the 
experiment.  The  nonlinearity  was  representing  the  relay  characteristic  on  the  on-off  valves. 

In  the  Random  Input  Describing  Function  [20]  analysis  one  assumes  that  the  input  u(t)  to  the 
nonlinear  element  is  of  the  form  u(t)  =  B  +  Asm(u>t  +  0)  -f-r(t),  where  B,  A ,  u)  are  unknown  constants, 
8  is  an  arbitrary  initial  phase,  and  r(t)  is  a  Gaussian  process  with  a  standard  deviation  a.  The  output 
of  the  nonlinear  element  y(t)  —  f(u(t))  is  approximated  as 


Va(t)  =  NbB  -I-  NAAsm(ujt  +  9)  +  NRr(t ), 


(1.3) 


where  the  individual  gains  (called  Describing  Functions  of  the  nonlinearity)  are  obtained  by  minimiza¬ 
tion  of  the  variance  of  the  residual: 


NB(B,A,v)  =  $E[f{u{  0))]  = 

J?  d9fZo  drf(B  +  Asin{9)  +  r)  exp(-^r), 


(2^)5 


oB 


Nr(B,  A,  a)  =  £E[f(u( 0))r(0)]  = 

/o2’  d0  f-oo  drf(B  +  AsinW  +  r)r  exP(~5^)> 


Na(B,  A,  a)  =  %E[f  (u(0))  sin(0)]  = 

Jo2’  d6f-oo  drf(B  +  AsMe)  +  r)  sin(0)  exP(-^r)- 


(1.4) 


(1.5) 


(1.6) 


([20],  p.  371).  For  the  no-noise  and  no-bias  case  where  r(t)  =  0  and  B  =  0,  the  last  formula  reduces 
to  the  standard  sinusoidal  input  describing  function  gain. 

This  framework  easily  allows  us  to  study  the  response  of  the  pressure  output  in  either  open  or 
closed  loop  thermoacoustic  system  as  an  output  of  a  stable  noise  driven  system  (when  A  =  0)  or 
as  self-excited  oscillation  (when  A  >  0)  or  some  combination  of  the  two.  We  assume  the  Gaussian 
component  r<(t)  of  the  input  disturbance  has  Power  Spectral  Density  $n{jw)  and  write  corresponding 
equations: 

1.  Stable  driven  system 

Go(0) 


B  = 


A  = 


$x  xCM  = 


1  +  Nb(B,A,<t)Go{0) 

Go  CM 


Bi 

Ai 


l+NA(B,A,a)G0(ju>) 

Go  CM  |2 


'l  +  NR(B,  A,o)Go(ju) 
1 

o1  =  —  J  $xx(juj)duj. 
2.  Self-excited  oscillations  with  driving  noise 

Go(0) 


*«CM 


B  = 


$x*CM  = 


1  +  Nb(B,A,ct)Go{0) 

1  +  Na(B,  A,<j)Go(ju)  —  0 
Go  CM 


Bi 


1  +  Nr(B,  i4,a)Go(M 

2  =  h  Loo  ^xx^duJ' 


!$«CM 


(1.7) 

(1.8) 

(1.9) 

(1.10) 

(1.11) 

(1.12) 

(1.13) 

(1.14) 
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Figure  1.2:  The  amplitude  of  the  limit  cycle  in  the  presence  of  the  Gaussian  noise 


Figure  1.3:  Describing  function  gains  in  the  presence  of  Gaussian  noise 

For  relay  nonlinearity,  the  amplitude  of  the  limit  cycle  can  be  found  from  solving  the  loop  equation 

(1.12) 

l  +  NA(A,a)G0{juj)  =  0.  (1.15) 

In  Mehta  et  a/.,  IFAC  2002,  we  prove  that  for  a  relay  nonlinearity: 

1.  For  a  — ¥  0,  the  appearance  of  the  limit  cycle  stabilizes  the  loop  with  respect  to  the  noise  thereby 
yielding  a  bounded  input-output  response  (for  the  Gaussian  noise  driver)  as  solution  of  equation 
(1.13)  and 

2.  for  a  — >  oo,  large  noise  stabilizes  the  loop  with  respect  to  the  limit  cycle  thereby  causing  the 
limit  cycle  to  disappear  and  the  system  to  behave  as  a  stable  noise  driven  system. 

More  precisely,  Figure  1.2  shows  that  the  presence  of  noise  (cr  >  0)  leads  to  a  reduction  in  the  amplitude 
of  this  limit  cycle  and  at  a  critical  positive  value  of  a  =  cr0,  the  limit  cycle  disappears  (A(ao)  =  0 
for  values  of  a  >  <to).  Figure  1.3  shows  the  gains  Nr(ct)  and  Na(&)  as  function  of  cr.  For  the  values 
of  a  where  limit  cycle  is  present,  Nr  monotonically  increases  between  0  and  <7o  and  decreases  for 
values  of  a  >  ctq.  We  also  showed  that  the  feedback  interconnection  of  G^joS)  and  Nr(A(o),(j)  is 
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linearly  stable  for  all  a  ^  ao  and  the  largest  loop  gain  occurs  at  the  critical  value  oq  where  the  loop  is 
arbitrarily  close  to  destabilization  (eigenvalues  on  the  imaginary  axis).  For  values  of  a  away  from  <7o, 
the  eigenvalues  move  in  to  the  LHP  thereby  ensuring  asymptotic  stability  for  all  a  <7o. 

Therefore,  except  for  a  critical  a  =  cr0,  the  Fourier  transform  of  the  Gaussian  component  of 
combustor  pressure  may  be  approximated  as 

PgO'w)  =  Go(jv)S(jw,  A, a)ni(jw),  (1.16) 


where  (jco)  is  the  Fourier  transform  of  the  driving  disturbance  and 


S{ju),A,a)  = 


1  +  Go(ju>)NR(Q,  A ,  a)Gc(jw) 


is  a  modified  sensitivity  function  that  depends  on  the  magnitude  of  the  limit  cycle  A  and  standard 
deviation  a  of  the  Gaussian  component  at  the  input  of  the  relay  nonlinear  element.  This  result  allows 
to  extend  the  results  of  the  standard  linear  fundamental  limitations  analysis  presented  in  the  papers 
[2,  3,  4]  to  the  case  of  control  with  on-off  actuators.  In  particular,  it  can  be  shown  that  the  supression 
of  the  pressure  oscillations  in  the  central  frequency  band  of  the  natural  combustor  resonance  (which 
requires  |5(ju;,  A,  a)  \  <  1)  will  lead  to  amplification  of  the  pressure  oscillations  in  an  adjacent  frequency 
band  (|5(jfu;,  A,  cr)|  >  1),  which  leads  to  peak-splitting.  The  norm  used  here  is  the  H0 0  norm. 

Analysis  presented  here  indicates  that  the  peaking  phenomenon  defined  as  excitation  of  oscillations 
with  closed-loop  control  is  to  a  large  extent  inevitable  for  combustion  processes  with  large  delay 
controlled  with  actuators  of  limited  bandwidth. 

In  Teerlinck  et  a/.,  2005,  we  applied  the  nonlinear  frequency- domain  framework  described  above 
to  explain  results  in  active  fuel  control  experiment  in  800kW  three-flameholder  rig.  The  experiment 
involved  on-off  valves  and  severe  peak-splitting  was  observed.  The  optimal  amount  of  fuel  w as  proposed 
to  reduce  the  peaking  and  to  provide  acceptable  attenutation  of  oscillations. 


1.2  Spectral  Balance 

A  frequency-domain  framework  for  analysis,  computations,  and  uncertainty  propagation  in  nonlinear 
systems  driven  by  broad-band  disturbances  was  introduced  and  illustrated  in  a  simple  example  of 
a  nonlinear  system  that  exhibits  noise-induced  transitions  between  two  locally  stable  equilibria.  An 
approximate  and  iterative  spectral  balance  (including  determination  of  equilibria)  is  solved.  The 
solution  of  the  approximate  spectral  balance  is  used  to  reformulate  the  original  model  using  a  loop 
transformation  so  that  an  iterative  procedure  for  finding  the  spectrum  of  the  output  converges  to  the 
true  spectrum  of  the  solution.  The  work  is  presented  in  paper  by  Banaszuk  and  Mehta,  CDC  2004. 

Many  industrial  flows  involve  complex  interactions  of  acoustic  waves,  vorticity,  fuel  transport,  and 
chemical  reactions.  The  control  objective  often  is  to  create  beneficial  non-equilibrium  dynamics  with 
control.  Examples  include  control  of  flow  separation  and  mixing  enhancement.  In  this  paper  we  intro¬ 
duce  a  frequency  domain  framework  for  analysis  and  non-equilibrium  control  design  for  a  large  class  of 
models  of  physical  phenomena  involving  multiple  oscillatory  modes  coupled  through  nonlinear  terms, 
transport  delay,  and  driven  by  broad-band  disturbances.  While  motivated  by  specific  problems  arising 
in  military  aeroengines,  the  methods  will  be  applicable  to  large  class  of  distributed  dynamical  systems 
involving  oscillatory  dynamics  with  nonlinear  cross-coupling,  saturated  nonlinearities,  transport  delay, 
and  broad-band  disturbances. 

The  spectral  balance  framework  that  we  propose  generalizes  the  standard  harmonic  balance  and 
Gaussian  signal  balance  in  feedback  systems  [31,  20].  The  framework  is  introduced  and  illustrated  in  an 
example  of  a  nonlinear  system  that  exhibits  noise-induced  transitions  between  two  stable  equilibria. 
The  example  presented  is  a  scalar  model  with  cubic  nonlinearity  after  pitchfork  bifurcation  driven 
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by  a  broad- band  disturbance.  An  approximate  and  iterative  spectral  balance  of  the  constant  and 
broad-band  signals  (including  determination  of  equilibria)  is  solved.  The  solution  of  this  approximate 
spectral  balance  is  used  to  reformulate  the  original  model  using  a  loop  transformation  so  that  an 
iterative  procedure  for  finding  the  spectrum  of  the  output  converges  to  the  true  spectrum. 

Consider  a  model  of  a  lightly  damped  stable  linear  system  with  transfer  function  Go(juj),  in  a 
feedback  loop  with  a  static  nonlinearity  /(•),  subject  to  a  driving  disturbance  n(£)  with  the  Fourier 
transform  N(ju) ).  An  uncertainty  in  the  model  is  represented  by  an  (in  general  nonlinear)  operator 
A(-)  in  a  feedback  loop  around  the  nominal  model  The  model  equations  are 


Nominal  Model 


Figure  1.4:  The  model  structure 

X(ju)  =  G0(ju)(N(ju)  -  X(jcj))  -  Y(jco))  (1.18) 

y(t)  =  f(x{t))  (1.19) 

where,  X(-)  =  Tx('),  Y(-)  =  ^y(-),  and  N(-)  =  JFn(*),  are  the  Fourier  transforms  of  the  corresponding 
temporal  signals.  We  assume  that  nonlinear  mapping  /(•)  is  Lipshitz  on  each  bounded  set.  The 
equation  (1.19)  can  be  represented  in  the  frequency  domain  as 

Y(juj)  =  f(X(ju>))  :=  JF/tJF-1  X(jw)).  (1.20) 

Now,  the  feedback  system  (1.18)— (1-19)  can  be  represented  as 

X{ju>)  =  G0(jv)(N(ju>)  -  f(X(jv)  -  A (MX(ju>)).  (1.21) 

Note  that  for  the  hnear  case  f{x)  =  0,  A (juj,X(juj))  =  A (ju))X(juj)  the  mapping  of  the  uncertainty 
A(ju;)  to  the  output  of  the  system  is  given  by  the  formula 

X(jw)  =  (I  +  Go(juj) A(jco))  lGo(ju))N(juf).  (1.22) 

involving  the  sensitivity  function  (7  +  Go(i^)A(jcj))_1.  Note  that  the  frequency  domain  representa¬ 
tion  greatly  simplifies  the  uncertainty  propagation  analysis. 

1.  Uncertainty  propagation.  The  sensitivity  function  (7  +  Go(j'cj)A(ja;))“1  allows  to  explicitly  map 
the  probability  distribution  of  the  uncertain  parameters  contributing  to  A (ju)  to  the  probability  dis¬ 
tribution  of  the  output  Y(jcj). 
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2.  The  frequency  domain  representation  greatly  accelerates  computation  of  this  mapping.  Note  that 
only  the  algebraic  calculations  need  to  be  performed  in  evaluating  the  formula  (1.22).  In  contrast, 
a  time  domain  counterpart  of  the  (1.22)  would  require  evaluation  of  the  convolution  integrals  over 
long  period  of  time.  Computation  of  large  linear  systems  with  significantly  separated  time  scales  is 
cumbersome  in  time  domain,  as  the  shortest  time  scales  detemine  the  time  step  size,  while  the  longest 
time  scales  determine  the  total  time  of  simulation.  Moreover,  in  the  time  domain  formulation  one 
needs  to  wait  for  transients  to  subside,  which  is  an  issue  when  dealing  with  lightly  damped  dynamics. 
There  are  additional  benefits  of  the  frequency  domain  representation  in* the  case  when  Gq(Jlj)  contains 
time  delays. 

3.  Tools  from  the  robust  linear  control  theory  allow  to  handle  dynamic  uncertainty  in  case  when 
only  the  bounds  on  the  uncertain  operator  are  known  [15]. 

Appart  from  application  to  the  uncertainty  propagation,  the  frequency  domain  formulation  allows 
to  study  fundamental  limitations  of  achievable  control  performance  using  methods  of  the  complex 
analysis. 

The  spectral  balance  approach  retains  the  advantages  of  the  linear  sensitivity  function  framework: 
explicit  formulas  mapping  uncertainty  to  the  output  and  the  speed  of  computation.  We  will  begin 
with  the  particular  case  of  system  in  Figure  1.4  using  the  fixed  point  formulation  (1.21).  To  introduce 
the  spectral  balance  framework  we  will  consider  the  case  without  uncertainty  shown  in  Figure  1.5  with 


Figure  1.5:  The  model  structure 

the  corresponding  fixed  point  formulation  of  the  spectral  balance  equation  given  by  (1.23) 

X(jco)  =  G0(ju)(N(jw)  -  f(X(juj)).  (1.23) 

Note  that  the  spectral  balance  framework  generalizes  the  standard  harmonic  balance  (where  the  input 
signal  n(-)  is  periodic,  or  when  the  dynamics  has  limit  cycles)  and  Gaussian  signal  balance  (where  the 
input  signal  n(*)  is  a  Gaussian  broad  band  signal)  in  feedback  systems  [31,  20]. 

We  assume  that  the  dynamics  of  (1.19)  is  globally  bounded  and  that  there  is  an  attractor.  Eventu¬ 
ally  we  intend  to  introduce  a  spectral  balance  framework  for  the  class  of  bounded  power  signals  on  an 
infinite  time  interval.  In  this  paper  we  restrict  the  attention  to  the  space  of  L<i  signals  on  the  interval 
[0,  T],  where  T  is  large  relative  to  the  slowest  time  scale  in  the  system.  The  induced  operator  norms 
are  the  Hqq  norms. 

A  sufficient  condition  for  existence  of  a  unique  solution  of  the  spectral  balance  equation  (1.23)  is 
\\GoUu)(f(X2(ju>))  -  /(Xitiw))) II  <  II X2(ju>)  -  Xi(ju,)||.  (1.24) 
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For  all  Xi(ju>)  in  L2[0,T].  Note  that  in  this  case  a  unique  solution  to  (1.23)  exists  (by  applying 
the  Banach  Contraction  Mapping  Theorem  [28]).  Moreover,  the  approximate  solution  of  the  spectral 
balance  equation  can  be  obtained  by  successive  approximations  using  the  formula 

Xi+l(ju)  =  G0{juj)(N(jw)  -  f(Xi(ju)).  (1-25) 

with  an  arbitrary  initial  condition. 

Note  that  if  the  condition  (1.24)  is  satisfied  for  all  Xx{ju>)  in  a  closed  set  B  in  L2[0,T]  that  is 
invariant  for  the  mapping  Go(jui)(N (juj)  —  /(•)),  one  can  approach  a  solution  of  the  spectral  balance 
equation  (1.23)  in  B  using  (1.25)  with  Xo(ju>)  £  B. 

1.2.1  Loop  transformation 

A  sufficient  condition  for  (1.24)  is  the  small  gain  condition  for  the  feedback  loop  in  (1.5).  However,  even 
if  the  condition  (1.24)  is  not  satisfied,  which  would  be  the  case  if  the  loop  gain  is  large,  one  can  attempt 
to  enforce  the  condition  (1.24)  for  an  equivalent  feedback  system  to  (1.5)  by  a  loop  transformation. 
An  example  of  a  linear  loop  transformation  is  shown  in  Figure  1.6.  Here  H{ju))  is  an  arbitrary  stable 


Unear  operator, 

Gi(ju)  :=  (I  +  H{ju;)Go(joj))-1Go(ju). 


(1.26) 


arid 

h(X{ju>))  ■■=  f(X{juj))  -  H{ju>)X{ju>).  (1-27) 

The  spectral  balance  condition  for  the  system  in  Figure  1.6  is 

X(ji j)  =  Gi{jco){N{ju)  -  fi{X{ju))).  (L28) 

Note  that  a  sufficient  condition  for  the  contraction  condition  for  the  transformed  spectral  balance 

(1.28)  is  . 

||Gi(^)(/i(X2(ju;))  -  A(X1(yw)))||  <  \\X2(ju)  -  (1.29) 

If  the  nonlinear  part  of  the  loop  in  Figure  1.5  has  a  stabilizing  effect,  the  role  of  the  operator  is 

to  reduce  the  ffoo  gain  of  the  nonlinear  part  of  the  loop  and  increase  the  contractivity  of  the  Unear  part 
of  the  loop.  More  precisely,  the  gain  of  the  nonlinear  operator  f(X(juj))  is  reduced  by  subtraction  of 
a  linear  approximation  of  f(X{juj))  and  the  approximate  Unear  operator  is  incorporated  in  the  linear 
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part  of  the  modified  loop.  Thus,  a  good  choice  of  H(juj)  is  the  one  that  minimizes  \\fi(X(juj))\\  for 
X(ju ;)  representing  the  solution  of  the  spectral  balance  equation  (1.28). 

Of  course,  the  minimization  of  \\fi(X(ju))\\  requires  knowledge  of  X(ju)  itself,  which  is  exactly 
the  solution  of  the  spectral  balance  equation  that  we  seek.  Note  that  we  are  interested  in  the  case  of 
system  (1.19)  having  non-equilibrium  attractors  or  subject  to  a  large  driving  disturbance,  so  that  an 
approximation  of  nonlinear  operator  f(X(ju)))  by  its  linearization  at  (X(juj)  =  0  is  not  suitable. 

The  key  idea  introduced  in  this  paper  is  to  proceed  in  the  following  three  steps: 

1.  Find  an  approximate  solution  Xappr(juj)  close  to  X(jco).  For  this,  the  describing  function  tech¬ 
niques,  both  for  harmonic  and  random  Gaussian  signals,  can  be  utilized.  In  fact,  as  we  will  show  in 
the  next  section,  it  may  be  not  necessary  to  find  an  approximate  solution  Xappr(j u>),  but  only  few 
parameters  describing  such  a  solution,  like  its  time  average  and  the  average  power. 

2.  Utilize  the  knowledge  of  Xappr(juj)  (or  the  parameters  describing  it)  to  find  a  linear  transforma¬ 
tion  H(ju )  that  minimizes  (or  at  least  reduces)  \\fi(Xappr(jcj))\\  =  1 1 / {Xappr (jw ))—H (ju))Xappr (jw) \ \ . 

3.  Use  H(jcj)  to  define  the  loop  transformation  (1.26)-(1.27).  If  the  contraction  condition  (1.29) 
is  satisfied  on  a  closed  set  B  in  L2[0,  T]  that  is  invariant  for  the  mapping  Gi{ju>)(N(ju ))  —  /i(-)),  one 
can  approach  a  solution  of  the  spectral  balance  equation  (1.23)  in  B  using  the  iterative  process  defined 
by 

Xi+l(ju>)  =  Gi(ju)(N(ju>)  -  MXiUuj))  (1.30) 

starting  with  an  arbitrary  Xq  (jw)  €  B . 

At  present  it  is  not  clear  under  what  general  conditions  the  procedure  described  above  will  result 
in  finding  solutions  to  the  spectral  balance  equations.  In  the  next  section  we  will  show  one  example 
of  a  system  with  nontrivial  dynamics,  for  which  that  the  procedure  yields  the  desired  result. 

1.2.2  Example 

Consider  the  equation 

x(t)  +  ax(t)  +  bx3(t)  =  n(t).  (1-31) 

Here  we  assume  that  a  <  0,  b  >  0,  and  the  input  signal  n(  )  has  zero  mean  and  flat  spectrum 
| N (j uj)  |  =  <ii  for  all  u;.  In  the  sequel  we  will  refer  to  the  input  signal  n(*)  as  noise,  even  though  we 
emphasize  that  in  this  paper  we  only  consider  the  deterministic  case.  Note  that  for  a  =  0  the  system 
(1.31)  undergoes  a  pitchfork  bifurcation  and  the  equilibrium  x  =  0  becomes  unstable  for  all  a  >  0. 
Two  locally  stable  equilibria  occur  at  x  =  and  x  =  —  yj—  Note  that  for  a  small  value  of  o\  the 

solution  x(t)  will  be  close  to  one  of  the  stable  equilibria.  For  some  higher  value  of  of  the  solution  x(t) 
will  be  transitioning  from  a  neighborhood  of  the  one  of  the  stable  equilibria  to  the  other,  as  shown  in 
Figure  1.7.  Note  that  the  spectral  balance  equation  for  (1.31)  is 

X(jw)  =  (N(ju)  -  f(X(ju>)),  (1.32) 

JUJ  +  a 

where  f(x)  :=  6x3.  Note  that,  since  a  <  0,  the  linear  operator  is  unstable,  and  the  contraction 
condition  (1.24)  is  not  satisfied.  However,  the  nonlinear  operator  f(x)  has  a  stabilizing  effect,  so  that 
we  can  attempt  to  transform  the  loop  to  an  equivalent  one,  for  which  the  contraction  condition  is 
satisfied,  as  described  in  Section  1.2.1. 

Let  x  :=  ^  Jq  x(t)dt  denote  the  time  average  of  x(t)  and  let  xf(t)  :=  x(t)  —  x  denote  the  deviation 
of  x(t)  from  its  average  value.  Let  o\  ^  x'(t)2dt  denote  the  mean  power  of  xf(t).  In  what  follows 
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Typical  solution 


Tima  trace,  a=-1 ,  b=1,  nowe  power=0. 99472 


Tuna, sac 


Figure  1.7:  Typical  solution  of  system  (1.31):  noise-induced  transitions  between  two  stable  equilibria 


we  will  compute  approximate  value  of  x  and  a2  by  solving  approximate  spectral  balance  equations.  By 
taking  the  time  average  of  (1.33)  and  using  the  fact  that  the  average  of  (x  +  x'(£))3  is  x3  +  3ba2  we 
obtain 

x(a  4-  3ba2  +  bx2)  =  ft  =  0.  (1.33) 

Subtracting  (1.33)  from  (1.31)  and  re-arranging  terms  yields 

x'(t)  +  (a  +  h)x{t)  +  fx(x'(t))  =  n(t),  (1.34) 


where 

h  :=  bcr 2  4-  3bx2 

fl{x'(t))  :=  b(x'2  -  o2x)(x'  +  35). 


(1.35) 

(1.36) 


To  find  approximate  values  of  x  and  a2  we  will  neglect  the  term  h(x,  axi  x'(i))  in  (1.34)  and  solve  the 
equation 


xr(t)  4-  (u  4~  bd2  *+■  3 bx2)xf  (t)  —  n(t ). 


(1.37) 


Note  that  for  fixed  values  of  x  and  a2  (1.37)  is  a  linear  equation  that  can  be  solved  in  the  frequency 
domain  as 


X'Uu) 


Njjco) 

ju>  +  a  +  ba2  -h  3 bx2' 


(1.38) 


For  a  moment  we  assume  that  the  values  of  x  and  a2  are  such  that  a  4-  ba2  4-  3bx2  >  0,  so  that 
the  transfer  function  ju+a+ba 2+3frx3  ^  This  assumption  will  be  verified  after  x  and  a2  are 

calculated. 

Now  we  obtain  the  closure  equation  for  a2  by  integrating  the  absolute  values  of  both  sides  of  (1.38) 
over  all  frequencies 


J_  f°°  * _ <H _ .2  . 

2n  J-oo  jv  4-  a  4-  ba2  4-  3 bx2 


(1.39) 
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The  equations  (1.33)  and  (1.39)  form  an  approximate  spectral  balance  for  the  system  (1.31.  The  integral 
in  (1.39)  can  be  analytically  evaluated  so  that  we  can  write  the  following  equation 


rr2  - 

(1.40) 

°x  ~~ 

Now,  solving  (1.33)  and  (1.40)  we  obtain 

2  (a  +  ba2  +  3  bx2) 

2  0,2 
ai  K  4 

(1.41) 

«5 = Ba  -  Vi  -  ¥). 

(1.42) 

ry 

(1.43) 

o  az 

CTi  >  46  ^ 

(1.44) 

(1.45) 

x  =  0. 

(1.46) 

Figure  1.8  graphically  represents  the  solution  to  the  approximate  spectral  balance  equations  as  function 
of  the  input  power  a*  for  a  =  —1,  6  =  1.  The  solutions  for  x  and  a2  have  a  natural  intepretation. 


Approximate  spectral  balance 
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Figure  1.8:  Solutions  to  approximate  spectral  balance  equations  as  function  of  the  input  power 
2 

For  tii  <  there  are  two  values  of  the  time  average  x  close  to  the  no-noise  equilibria  that  can  be 
attained.  The  value  of  ax  (that  could  be  interpreted  as  standard  deviation  of  x(t))  is  small,  so  that 
the  solution  x(t)  stays  close  to  an  equilibrium  solution  and  does  not  transition  to  the  neighborhood  of 
the  other  equilibrium.  Above  the  critical  value  of  the  input  power  ai  =  ^  the  solutions  x(t)  deviate 
from  the  stable  equilibria  far  enough  to  transition  between  the  neighborhoods  of  the  both  equilibria. 
Since  the  transitions  back  and  forth  can  occur,  x  =  0  becomes  the  mean  and  the  standard  deviation 
ax  is  close  to  the  distance  from  the  new  mean  x  =  0  to  the  value  where  the  solution  x(t)  spends  most 
of  the  time:  close  to  the  no-noise  equilibria  of  (1.31). 
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We  will  now  use  the  values  of  x  and  a\  given  by  (1.8)  to  perform  a  loop  transformation  as  described 
in  Section  1.2.1.  More  precisely,  we  will  solve  the  perturbation  equation  (1.34)  in  the  frequency  domain 
using  the  fixed  point  formulation 

X'(juj)  =  -±-—(N(jw)  -  A (X'(ju)).  (1.47) 

jw  +  a  4-  h, 

It  can  be  easily  verified  that  a  4-  h  >  0  for  x  and  given  by  (1.46).  Analytic  verification  of  the 
contraction  condition  for  the  operator  ,^+a+A (N(juj)  —  /i (•)  is  difficult.  Therefore  we  will  assume  that 
the  contraction  condition  is  satisfied  and  proceed  with  an  iterative  solution  to  (1.47)  using 

X’i+i(jw)  =  -±—-(N(jw)  -  A (Xiticj))  (1.48) 

ju)  -b  a  -t -  ti 

with  X*q  (ju)  =  0  and  verify  the  contraction  condition  numerically.  To  illustrate  and  verify  this 

2 

procedure,  a  numerical  solution  of  (1.31)  for  a  =  —1,  6  =  1,  and  a  >  was  obtained.  The  spectrum 
N(juj)  of  the  noise  from  the  time  domain  simulation  was  saved  and  used  in  the  formula  (1.48). 
Figure  1.9  shows  an  excellent  agreement  of  the  spectrum  X^juj)  from  the  time  domain  simulation 
and  the  spectrum  X'io(juj)  from  the  iterative  procedure  (1.48)  after  10  iterations.  Figure  1.10  shows 

Simulation  and  approximation 
after  10  iterations:  spectra  (jft) 

Spectra,  a=-1 ,  b=1,  noise  power=3.9789 


Frequency,  Hz 

Figure  1.9:  Solution  to  iterative  spectral  balance  equations:  spectra 

comparison  of  the  time  traces  of  the  solutions  of  (1.31)  obtained  by  the  time  domain  simulation  and 
by  the  iterative  spectral  balance  and  the  inverse  Fourier  transform.  Figure  1.11  shows  decay  of  the 
power  of  the  approximation  error  Xf(ju)  —  X'i(ju))  normalized  by  the  power  of  X'(juj)  as  a  function 
of  iteration  step  i.  Finally,  Figure  1.12  shows  the  contraction  rate  1 ^ ^  35  a  func!i°n 
of  iteration  step  i.  This  verifies  the  contraction  at  the  rate  of  about  0.8  was  indeed  achieved  by  the 
loop  transofrmation  involving  solution  of  the  approximate  spectral  balance. 

In  the  cases  that  were  examined  obtaining  an  approximate  solution  of  (1.31)  using  the  formula 
(1.48)  was  orders  of  magnitude  faster  that  the  time  domain  simulations  using  Simulink. 
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Simulation  and  approximation 
after  10  iterations:  time  traces 

Time  traces,  a—1 ,  b=1 ,  noise  power=3.9789 


- Time  simulation 

-  Iterative  Spectral  Balance 

—  Equilibria 
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Figure  1.10:  Solution  to  iterative  spectral  balance  equations:  time  traces 

1.2.3  Conclusion 

A  frequency-domain  framework  for  analysis,  computations,  and  uncertainty  propagation  in  nonlinear 
systems  driven  by  broad-band  disturbances  was  introduced  and  illustrated  in  a  simple  example  of  a 
nonlinear  system  that  exhibits  noise-induced  transitions  between  two  stable  equilibria.  The  spectral 
balance  framework  generalizes  the  standard  harmonic  and  Gaussian  signal  balance  in  feedback  sys¬ 
tems.  The  application  example  presented  is  a  scalar  model  with  cubic  nonlinearity  after  pitchfork 
bifurcation  driven  by  a  broad-band  disturbance.  An  approximate  and  iterative  spectral  balance  (in¬ 
cluding  determination  of  equilibria)  is  solved.  The  solution  of  the  approximate  spectral  balance  is 
used  to  reformulate  the  original  model  using  a  loop  transformation  so  that  an  iterative  procedure  for 
finding  the  spectrum  of  the  output  converges  to  the  true  spectrum  of  the  solution.  The  future  work 
will  involve  more  carefull  study  of  the  function  spaces  suitable  for  the  spectral  balance  formulation 
and  obtaining  some  analytic  sufficient  conditions  for  the  contraction. 

1.3  Adaptive  Control  of  Flow  Phenomena  in  Aeroengines 

In  paper  by  Banaszuk  et  al ,  Automatica  2004,  we  described  adaptive  control  scheme  for  control  of 
oscillations  with  unknown  frequency  and  amplitude  and  its  application  to  control  of  thermoacoustic 
instabilities.  The  original  submission  was  supported  by  the  previous  AFOSR  grant,  but  the  subsequent 
revisions  were  supported  by  the  current  grant. 

In  paper  by  Krstic  and  Banaszuk,  Control  Engineering  Practice  2003,  we  considered  the  problem 
of  stabilization  of  a  class  of  MIMO  LTI  systems  arising  in  models  of  various  instabilities  in  jet  engines. 
The  problem  was  motivated  by  control  of  flutter,  stall,  and  thermoacoustics  in  military  engines.  These 
instabilities  often  manifest  themselves  as  oscillations,  contaminated  by  noise.  They  are  often  caused 
by  coupling  of  several  resonant  modes  (structural,  acoustic,  of  vortical)  with  time  delays  present  in 
the  physical  process  that  couple  the  resonant  modes.  Often  the  control  input  is  also  subject  to  delay. 
Possible  applications  of  the  results  in  the  paper  include  control  of  compressor  blade  flutter,  rotating 
stall,  and  aeroacoustic  instabilities  (coupling  of  acoustic  waves  with  vortex  shedding  from  stator  vanes). 
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Figure  1.11:  Solution  to  iterative  spectral  balance  equations:  relative  approximation  error  for  210  FFT 
points 

Uncertain  parameters  abound  in  these  problems:  unknown  or  varying  natural  frequencies,  un¬ 
certain  delays  due  to  poorly  understood  physical  phenomena  governing  these  processes,  uncertain 
coupling  between  modes  of  oscillation,  and  of  course,  uncertain  high  frequency  gains  and  delays  of 
actuators.  In  this  paper  we  approach  a  class  of  such  models  using  the  tools  of  adaptive  control. 

Consider  the  model  of  the  form 

yi  +  £ii£i  +  miyi  +  Cn2/i(t  —  Tli)  +  Ci22/2(t  —  T12)  =  9n^2{t  -  rCy\)  +  hnx i  (1.49) 
92  -I-  C222/2  +  r?222/2  +  C22I/2  —  T22)  4-  C212/I (*  “  T21)  =  922^1  (t  —  TCy 2)  +  h>22X2  (1.50) 

where  y\  and  xj2  are  temporal  coefficiants  of  the  resonant  modes,  xi  and  X2  are  the  disturbance  inputs, 
and  the  parameters  & j ,  Vij  >  Cij  >  rij >  9ij ?  hij ,  rCy{  are  uncertain.  Such  a  model  is  common  in  case  where 
two  resonant  modes  with  close  resonant  frequencies  couple  though  a  physical  process  that  involves 
transport  delays.  For  example,  in  compressor  blade  flutter  the  variables  y\  and  y2  could  represent 
temporal  coefficients  of  blade  displacement  in  a  rotating  frame.  As  blades  move,  they  perturb  the 
flow.  In  turn,  the  flow  perturbations  affect  (with  some  delay)  the  blade  motion.  Because  the  blades 
have  airfoil  shape  and  the  mean  flow  has  swirl,  the  flow  reponse  is  not  axisymmetric  and  hence  can 
couple  the  resonant  modes.  The  cross  coupling  in  the  model  is  represented  by  the  terms  Cijyi{t  —  r^) 
with  i  /  j.  The  terms  CuViit  ~  ra)  represent  effect  of  flow  response  on  the  i-ih  mode,  which  can  be 
either  stabilizing  or  destabilizing. 

In  this  paper  we  were  interested  in  a  particular  case  of  strong  cross-coupling  of  identical  lightly 
damped  resonant  modes  represented  by  the  equations 

rn  +  VVi  +  CV2(t  -  t)  =  gu2(t-r)  +  hxi  (1.51) 

m  +  W2  -Cyi(f-  t)  =  -gui(t-r)  +  hx2-  (1.52) 

Such  an  interconnection  results  in  coupling  of  the  modes  yi  and  1/2-  K  the  latter  represent  temporal 
coefficients  of  resonant  standing  modes,  the  coupled  dynamics  often  represents  traveling  waves.  Note 
that  the  models  of  flutter  (Banaszuk  et  a/.,  IFAC  2002)  and  thermoacoustic  instabilities  on  annular 
domain  (Banaszuk  et  a/.,  CDC  2003)  presented  in  the  later  sections  of  this  report  are  of  a  similar 
form. 

The  adaptive  control  demonstration  in  such  model  is  described  in  details  in  Krstic  and  Banaszuk, 
Control  Engineering  Practice  2003. 
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Figure  L12:  Solution  to  iterative  spectral  balance  equations:  contraction  rate 

1.4  Uncertainty  Propagation  in  Complex,  Nonlinear  Interconnected 
Dynamical  Systems 

Under  DARPA  funding  (AFOSR  Contract  F49620-03-C-0035)  we  have  shown  that  the  uncertainty 
propagation  in  complex,  interconnected  dynamical  systems  can  be  performed  more  efficiently  by  de¬ 
composing  the  network  based  on  the  hierarchy  and/or  the  strength  of  coupling.  The  results  of  this 
research  are  summarized  in  CDC  2004  papers  by  Yarigonda  et  al. and  Huzmezan  and  Kalmar-Nagy 
presented  in  two  invited  sessions  on  Uncertainty  Propagation.  Some  basic  research  aspects  of  this  work 
were  analysed  in  more  detail  under  the  current  AFOSR  contract.  In  particular,  in  Varigonda,  CDC 
2004,  we  proposed  an  iterative  method  for  static  feedback  systems  to  obtain  the  probability  density 
of  the  output  from  that  of  the  input.  We  proved  the  convergence  of  the  proposed  method  under  the 
assumption  that  the  loop  operator  is  contractive.  The  method  was  illustrated  with  an  example.  It 
was  shown,  based  on  the  results  from  the  theory  iterated  random  functions,  that  the  method  extends 
to  the  case  when  additional  parametric  uncertainty  is  present  within  the  loop. 


19 


Chapter  2 

Control  of  Flow  Phenomena  in 
Military  Aeroengines 


2.1  Control  of  Mixing  in  Shear  Flows 

2.1.1  Control  of  Mixing 

Papers  by  Tadmor  and  Banaszuk,  IEEE  TCST  2002,  Wang  et  a/.,  Physics  of  Fluids  2003,  and  Noack 
et  at,  Physics  of  Fluids  2004  present  the  results  of  synergistic  use  of  Control  Theory  and  Dynamical 
Systems  methods  to  create  beneficial,  non-equilibrium  dynamics  in  low  dimensional  fluid  models.  The 
common  idea  is  to  use  control  to  enforce  a  periodic  behavior  in  fluid  velocity  that  creates  a  chaotic 
advection  field  for  the  fuid  particles. 

2.1.2  Control  of  Diffuser  Flow  Separation 

In  paper  by  Banaszuk  et  at,  AIAA  Reno  2003,  we  described  an  application  of  extremum-seeking  to 
adaptive  flow  control  in  a  subsonic  diffuser.  Specifically,  we  presented  results  of  an  experimental  study 
of  on-line  optimization  of  the  pressure  recovery. 

Separation  phenomena  occur  in  many  industrial  and  military  applications  including  external  flows 
such  as  flow  past  high  angle  of  attack  airfoils,  and  internal  flows  such  as  aggressively  expanding 
diffusers.  Consequently,  its  control  for  performance  improvement  has  received  widespread  attention. 
Various  means  for  delaying  the  onset  of  separation  have  been  proposed,  including  passive  and  active 
methods  [19].  The  use  of  periodic  oscillations  to  delay/reduce  the  extent  of  separation  in  airfoils 
was  investigated  (e.g.  see  [50]  ),  demonstrating  the  effectiveness  of  unsteady  blowing  in  controlling 
flow  separation.  Multi-frequency  open  loop  forcing  was  shown  to  create  and  enhance  interactions 
of  multiple  flow  structures  in  simple  free  shear  layers  [25]  Recently,  two-frequency  forcing  using  a 
synthetic  jet  actuator  was  shown  to  be  an  effective  way  of  increasing  diffuser  pressure  recovery  [40]. 
However,  it  is  difficult  to  predict  an  optimal  set  of  parameters  that  include  the  number  of  frequencies, 
relative  amplitude  and  phase  difference  between  the  forcing  frequencies  for  enhancing  performance, 
due  to  the  lack  of  an  analytical  or  modeling  method.  In  particular,  in  [40]  the  parameters  for  a  two- 
frequency  forcing  control  law  that  optimized  the  pressure  recovery  in  a  two-dimensional  diffuser  were 
found  manually. 

In  this  paper  we  present  a  method  for  automatic  tuning  of  parameters  of  a  multi-frequency  forcing 
flow  control  law  to  optimize  pressure  recovery  in  a  diffuser.  The  method  is  known  in  optimization 
and  control  theory  literature  as  extremum  seeking  [5]  In  extremum-seeking  control  one  adapts  the 
control  parameters  using  on-line  estimation  of  gradients  of  the  performance  metric  with  respect  to  the 
control  paramaters  by  introducing  small  probing  signals  on  top  of  (typically  slowly  varying)  control 
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parameters.  The  method  has  been  widely  applied  in  industry  [49,  1]  for  ” model-free”  optimization  of 
steady-state  of  many  industrial  processes. 

Fore  more  details  we  refer  to  Banaszuk  et  al. ,  AIAA  Reno  2003. 

2.1.3  Control  of  Pattern  Factor 

Key  performance  metrics  for  military  aeroengines  includes  pattern  factor,  controlled  primarily  by  the 
combustor  stochoimetry  and  the  degree  of  fuel-air  mixing.  In  a  typical  non-premixed  combustor, 
both  fuel  and  combustion  air  are  introduced  longitudinally  at  the  dump  plane,  and  swirl  is  generally 
utilized  to  mix  the  fuel  and  air  streams  together.  In  certain  combustor  designs,  additional  primary 
and  dilution  air  are  introduced  radially  through  circumferential  holes  located  along  the  combustor 
shell.  These  air-jets  in  cross  flow  not  only  provide  the  air  needed  to  control  the  stochiometry,  but 
also  generate  enhanced  fuel-air  mixing.  Therefore  the  proper  design  and  utilization  of  these  air-jets 
can  provide  a  means  toward  controlling  the  fuel-air  mixing,  and  enhancing  the  performance  metrics. 
Increased  mixedness,  in  particular,  can  provide  lower  pattern  factor. 

Control  of  jets  in  cross  flow  was  investigated  using  hierarchy  of  models  (including  high  fidelity 
CFD,  medium,  and  low  order  models)  and  described  in  Blossey  et  al .,  IUTAM  2001.  In  particular, 
low  order  mo  del- based  analysis  indicated  benefits  of  low  frequency  forcing  for  improved  mixing.  In 
the  current  funding  cycle  we  confirmed  the  model  prediction  in  experiment.  This  work  is  descibed  in 
paper  by  Narayanan  et  al .,  AIAA  Journal  2003. 

Furthermore,  the  benefits  of  jet  in  cross-flow  control  for  pattern  factor  reduction  were  demonstrated 
in  experiment  by  our  academic  partners  from  Luisiana  State  University.  The  experiments  were  not 
funded  by  AFOSR.  However,  the  UTRC  personnel  participation  in  this  joint  project  was  funded  by 
the  current  AFOSR  grant.  This  work  is  descibed  in  detail  in  Tuncer  et  a/.,  ASDME  2003.  Here  we 
provide  a  summary  of  this  work. 

The  effect  of  a  forced  dilution  air  jet  introduced  through  the  combustor  shell,  on  the  air-fuel  mixing 
in  the  combustion  chamber  has  been  investigated.  Thermocouple  based  temperature  measurements 
have  been  made  at  a  number  of  forcing  frequencies  in  the  range  of  lOO-HOOHz  and  blowing  ratios  in  the 
range  of  10-15.  The  open-loop  integral  flame  response  to  forcing  has  also  been  acquired  by  recording 
pressure  and  heat  release  spectra.  A  CH-radical  imaging  technique  is  used  to  provide  spatially-  and 
temporally-  resolved  information  about  the  heat  release  behavior.  The  results  exhibit  that  the  mean 
temperature  field  inside  the  main  reaction  zone  can  significantly  be  altered  as  a  consequence  of  air  jet 
modulation.  The  most  significant  effects  are  observed  by  forcing  at  vertical  locations  that  are  close 
to  the  dump  plane.  Enhancements  in  temperature  of  the  order  of  100-200  degrees  C,  and  reduction 
in  pattern  factor  of  the  order  of  10%  (e.g.,  from  1.13  to  1.03)  were  observed,  with  the  lowest  pattern 
factors  achieved  at  the  lowest  forcing  frequency  of  100  Hz. 

2.2  Modeling,  Analysis,  and  Control  of  Thermoacoustic  Instabilities 
in  Military  Engines 

In  the  past,  we  have  considered  reduced  order  models  of  combustion  instability  based  on  lumped 
parameter  modeling  as  described  in  [43,  26].  The  aim  of  this  Section  is  to  summarize  some  of  the 
more  recent  (and  ongoing)  research  whose  purpose  is  to  obtain  and  analyze  distributed  models  of 
thermoacoustic  instabilities  in  bluffbody  flameholder  annular  combustors.  These  instabilities  arise  on 
account  of  complicated  interactions  between  combustion  (flame  and  fuel)  dynamics,  vortex  dynamics 
and  acoustics  (see  Figure  2.1  for  a  schematic).  Our  strategy  thus  far  has  been  to  model  individual 
pieces  of  this  inter-connected  system  with  a  view  of  understanding  them  as  a  first  step  towards  gaining 
understanding  of  the  whole.  Before  summarizing  the  results,  we  provide  a  brief  summary  of  the 
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Figure  2.1:  Schematic  showing  interconnection  of  combustion  (flame-fuel),  vortex  and  fuel  dynamics 
in  thermoacoustic  problem. 

modeling  activity.  There  are  two  areas  of  modeling:  Thermoacoustic  models,  with  a  view  to 
understand  the  dynamics  of  coupled  system  and  Heat  release  models,  with  a  view  to  understand 
the  combustion  and  vortex  dynamics  (sans  acoustics).  For  the  thermoacoustic  modeling,  we  have  three 
versions  of  models  which  are  being  studied  for  the  annular  combustor  problem: 

3D  Linear  model  is  described  in  the  Section  2.2.1.  This  model  is  linear  and  assumes  that  dynamics 
arising  due  to  vorticity  and  flame  motion  are  neglected.  The  model  allows  us  to  isolate  the  role 
of  fuel  dynamics  (shown  to  be  equivalent  to  a  distributed  delay  for  distributed  but  fixed  flame)  in 
analyzing  thermoacoustic  instabilities.  The  linearity  hypothesis  allows  us  to  apply  linear  (fuel) 
control  methods  for  controlling  rotating  wave  combustion  instabilities  in  3D  annular  combustors. 
Additional  details  are  provided  in  Section  2.2.1  and  in  the  paper  [32]. 

2D  nonlinear  model  is  described  in  the  Section  2.2.2.  In  this  model,  a  3D  non- vortical  model  of 
thermoacoustic  instability  is  averaged  in  the  radial  direction  to  obtain  a  2D  model  for  multiple 
flameholders.  The  nonlinearities  as  well  as  flame  dynamics  are  retained.  The  non  vortical  model 
is  summarized  in  [33]  and  is  being  used  for  uncertainty  analysis  and  design  of  so-called  liners 
for  suppressing  combustion  instabilities.  In  Section  2.2.2,  a  version  of  this  model  which  includes 
vorticity  by  formal  superposition  of  2D  vortex  dynamics  is  summarized. 

Single  flameholder  model  is  the  same  as  the  above  2D  nonlinear  model  (including  vorticity  and 
flame  dynamics)  but  the  problem  is  simplified  to  study  only  a  single  flameholder  configuration. 
The  details  of  this  model  are  also  summarized  in  Section  2.2.2. 

In  addition  to  thermoacoustic  models,  we  have  also  concentrated  on  studying  the  heat  release  piece 
of  the  thermoacoustic  model  separately  (where  acoustics  is  neglected).  There  are  two  reasons  for 
studying  the  heat  release  piece.  One,  the  complexity  of  the  thermoacoustic  problem  resides  in  the 
heat  release  submodel  (where  complicated  interactions  occur  between  flame,  fuel  and  vorticity)  and 
two,  fuel  and  flow  control  aimed  at  modifying  heat  release  distribution  with  a  view  of  controlling 
combustion  instabilities  can  be  studied  effectively  with  these  models.  Instead  of  presenting  the  heat 
release  submodel  separately,  we  describe  three  studies  that  have  been  undertaken  with  the  purpose  of 
analyzing  and  controlling  heat  release  models.  In  Section  2.2.3,  we  present  a  vortex  model  developed 
to  study  the  physics  of  reacting  bluffbody  wake  dynamics.  The  model  as  presented  concentrates  on 
the  interaction  of  vorticity  and  flame  dynamics.  In  the  subsequent  two  Sections,  we  use  reduced  order 
modeling  approaches  to  better  understand  the  reacting  flow  dynamics  and  to  model  other  pieces  of 
the  heat  release  model.  In  particular,  in  Section  2.2.4  we  present  a  reduced  order  model  study  of 
flame- fuel  interaction  and  in  Section  2.2.5  a  study  of  vorticity-flame  interaction  (the  latter  study  is 
carried  out  to  explain  some  of  the  results  of  Section  2.2.3). 
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2.2.1  Three  dimensional  linear  model 

Thermoacoustic  instabilities  in  gas  turbine  and  rocket  engines  develop  when  acoustic  waves  in  com¬ 
bustors  couple  with  an  unsteady  heat  release  field  in  a  positive  feedback  loop.  We  consider  an  annular 
combustor  that  includes  a  circumferential  array  of  bluff  body  flame  holders  [48]  [16].  Flameholders 
extend  radially  from  inner  to  outer  diameter  of  the  annular  combustor.  A  cut  along  a  constant  radius 
surface  is  shown  in  Figure  2.2. 


Figure  2.2:  Bluff  body  flameholder  array,  fuel  source  surface  upstream  of  flameholders,  and  flame 
surface  downstream  of  flameholders. 


For  the  purpose  of  modeling,  we  assume  that  the  fuel  mass  fraction  defined  at  the  fuel  injection 
surface  xq (y,  z)  is  advected  downstream  to  the  fixed  but  distributed  flame  surface  x  =  g/i(y,  z)  by  the 
sum  of  the  mean  and  acoustic  perturbation  velocity  (without  diffusion).  The  mean  fuel  mass  fraction 
at  the  fuel  injection  surface  is  Y /(zo,y>  z)  =  >  where  Xf  ~  PfUf,  Xa  =  paUa  denote  the  flux 

and  pfj  pa  are  the  fuel  and  air  densities  and  & y,  Ua  are  the  velocities.  The  perturbation  fuel  mass 
fraction  (in  the  presence  of  acoustics)  is 


»/(*».».*.*>  =  7i(^)  (fss?  -  fgss?) 


(2.1) 


The  fuel- air  mixture  convects  to  the  fixed  flame  surface  x  —  gji(y,z )  and  the  heat  release  density  at 
the  flame  surface  is  obtained  as 


Q{x,y,z,t)  =  Fhr{Yf(x,y,z,t))'Yflame{x-gfi{y,z)))  (2.2) 

where  7 fiame{‘)  is  the  axial  heat  release  distribution  function  representing  the  flame  thickness,  and 
i?)lr(*)  describes  local  heat  release  as  function  of  local  fuel  mass  fraction. 

In  order  to  obtain  the  thermoacoustic  model,  we  define  relative  perturbations  of  pressure  and 
heat  release  as  p(x,y,z,t)  and  q(x,y,z,t)  :=  7~1  ,  where  7  is  the  ratio  of  specific 

heats.  We  also  assume  that  the  acoustic  velocity  perturbation  is  purely  potential,  i.e.,  u '(x,y,z,t)  = 
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V<f>(x,yyz,t)  for  some  smooth  scalar  <j){x,y,z,t)  called  the  velocity  potential  In  these  co-ordinates,  we 


obtain  a  linear  distributed  thermoacoustic  iiodel  (see  [32]  for  details)  as 

q 

0jp{x,  V >  z,  t)  +  u(x,  y ,  z)  ■  Vp( x,  y,  z,  t)  +  A <p{x,  y,  z,  t)  =  q(x,  y,  z,  t)  (2.3) 

Q 

— y,z,t)+  u(x,  y,  z)  ■  y,  z ,  t)  +  a2p(x , y,  z,  t)  =  r?(x,  y ,  z ,  t)  (2.4) 

Q 

^/(z,  y,  z,  t)  +  u(x,  y,  z)  •  Vy/(x,  y,  z,  t)  +  VY/(x,  y,  z)  •  V^(x,  y,  z,t)  =  0  (2.5) 

q{x,y,z,t)  =  F^Y f{x,y,z))'yflame(x  -  gfi(y,z))yf{x,y,z,t),  (2.6) 


where  driving  disturbance  (broad-band  noise)  rj(x ,  y,  z,  t)  represents  the  effect  of  local  turbulence.  The 
fuel  mass  fraction  boundary  condition  is  defined  on  the  fuel  injection  surface  as 

y/Qco(y>  z),3/,z,*)  _  ( ^(xo(y,z),p,z,t) _ 1 _ a<ft(x0(y,  z),y,z,f) 

Yf(x0(y,z),y,z)  \  uf{x0(y,z),y,z)  u (x0(y,z),y,z)  dx 

where  the  fuel  velocity  v!j  (xo(y,  z),  y,  z,  t)  is  the  control  variable.  The  first  term  on  the  right  hand  side 
of  (2.7)  represents  the  effect  of  fuel  control  action  and  the  second  term  represents  the  effect  of  acoustic 
velocity  perturbation. 

The  acoustic  boundary  conditions  are  provided  on  the  combustor  boundary  surface  in  terms  of  the 
normal  velocity  u^(x,  t)  =  V<£(x,  t)  *  n(x)  (where  fi(x)  is  the  normal  vector  to  the  boundary).  The 
acoustic  boundary  condition  serves  as  another  possible  control  input.  We  assume  that  the  acoustic 
boundary  conditions  are  described  by  a  local  admittance  relation  (described  here  in  the  frequency 
domain) 

K  (x,jw)  =  Gbc(x,ju)P(x,jw),  (2.8) 

(see  e.g.  [39])  for  xG<S,  where  S  denotes  the  boundary  surface. 

Reduced  order  model  for  control  of  thermoacoustic  instability  on  annular  domains 
The  fuel-air  mixture  is  responsible  for  the  burning  at  the  flame  and  the  subsequent  heat  release.  This 
heat  release  at  the  flame  surface  excites  the  acoustic  waves  in  the  combustor  volume.  The  acous¬ 
tic  waves  in  turn  travel  upstream  and  perturb  the  transport  of  the  fuel/air  mixture.  This  feedback 
coupling  can  lead  to  instability  if  the  driving  resulting  from  this  feedback  mechanism  dominates  the 
damping  resulting  from  absorption  of  the  acoustic  energy  at  the  boundary.  Control  over  fuel  rate  at 
the  fuel  injection  surface  xq (y,z),  control  of  the  shear  layer  dynamics  using  flow  control  at  the  flame- 
holders,  or  control  of  the  air  injection  at  the  combustor  boundary  can  provide  ways  of  influencing  the 
process  and  eliminating  instability.  The  control  could  be  provided  at  various  temporal  and  spatial 
scales. 

Now  we  introduce  a  reduced  order  model  (that  is  suitable  for  control  design)  obtained  from  the  ther¬ 
moacoustic  instability  model  presented  above.  The  model  is  also  suitable  for  optimization  of  the  control 
architecture.  In  order  to  obtain  model  reduction,  we  expand  the  pressure  and  potential  perturbations 
in  terms  of  the  acoustic  modes  {nfc(x)}fc=ii2,...  as  p(x,f)  =  £fcp*(f)IIjt(x),  <j>{x.,t)  =  £*  <^(t)nfc(x). 
and  apply  standard  Galerkin  procedure  involving  integration  by  parts  and  using  the  admittance  con¬ 
dition  (2.8)  (see  [32]  for  details)  to  obtain  a  two  mode  model  represented  in  the  frequency  domain 
as 

$i(jw)  0  0  —a2  0  Ni(jw) 

$2  O'w)  0  0  0  -a2  $2  (jw)  N2{ju>) 

JUJ  Pi  OH  Vo  o  o  Pi  CM  +  Qi(ju)-Vi(juj) 

.  P2(ju)  J  0  A2  0  0  J  [  p2(jw)  J  [  Q2(jw)  -  V2{jJ) 
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where  Qm(ju) )>  Nm{juj),  and  Vm{jui)^  denote  the  Fourier  transforms  of 


qm(t)  =  /  nm(x)q(x,  t)dx, 

Jv 

7]m(t)  =  /  nm(x)7?(x,t)dx, 

Jv 

=  /  nm(x)Un(x,t)dx, 

*/«S 


respectively,  and  Am  :=  ^  an^  ^  denotes  the  combustor  volume  and  boundary  sur¬ 

face  respectively).  The  fuel  velocity  ^(aroG/,  2),  y}  z,  t)  at  the  fuel  injection  surface  is  the  control 
variable.  We  assume  that  the  control  is  realized  using  N{nj  fuel  injectors  with  i  —  th  fuel  injector  pro¬ 
viding  fuel  mass  flux  equal  with  spatial  distribution  kj^{y,z)  (representing  initial  fuel  spread 

in  the  direction  perpendicular  to  the  mean  flow).  Thus,  we  will  represent  the  distributed  velocity  as 
W/OeoG/jz),?/,;*,^)  =  kf9i(y,z)wfti(t),  with  wj^t)  representing  the  control  inputs. 

The  closure  equations  to  (2.9)  are  given  by 

Vm(jw)  =  Gb£{ju))Pm{ju)  (2.13) 


-  ■  ITIJ 

Qm(ju)  =  GfZk(ju)<f>m(ju)  +  J2  GuJ9(jw)Wf,i(jw). 

k= 1  i- 1 

The  transfer  functions  in  the  above  expression  have  the  form 

Gm(jw)  :=  /5G6c(x,jw)|nm(x)|2rfx 

GrnfcO'w)  =  fw  ^^nm(9fi{y,z),y,z)F'(Yf{y,z))VY f{y,z)VUk(x,y,z)e~3 


_  _  .  9fi(y,x)-x 

jv)  =  fw  W^U^9fi{y,z),y,z)F'(Yf(y,z))VY f(y,z)vnk(x,y,z)e  3“  =<*.*>  dxdydg.  16) 
Grni  (jv)  =  fx0  uf(xolyyz)l,z)Urn(9fl(y, z), V, z)F'{Yf{y, z))e  3U1  ^  dydz  (2.17) 


While  the  reduced  order  frequency  domain  model  looks  deceivingly  simple,  it  is  in  fact  a  complicated 
infinite  dimensional  model,  as  the  heat  release  response  transfer  functions  include  a  distributed  delay. 

We  assume  that  pressure  measurements  p(xj,£)  =  (t)n*(xi)  at  one  or  several  locations  x* 

are  available.  We  also  assume  that  disturbance  terms  N\(ju ;),  Nzijw)  are  broad  band  uncorrelated 
stochastic  processes.  The  objective  of  the  feedback  control  is  to  reduce  the  gain  (Woo  or  W2)  from 
the  disturbance  terms  Ni,  N2  to  the  pressure  terms  Pi,  P2  to  guarantee  that  the  pressure  level 
is  below  acceptable  level.  Once  the  mean  fuel  mass  fraction  distribution  Yf(yyz)  and  control  fuel 
injection  (represented  by  choice  of  influence  functions  kf^y^z))  are  defined,  linear  control  laws  are 
straightforward  to  obtain.  However,  the  real  challenge  is  optimization  of  the  control  architecture. 
Namely,  one  would  like  to  select  the  mean  fuel  mass  fraction  distribution  Yf(y,z)  and  control  fuel 
injection  influence  functions  kfj(y,  z)  that  guarantees  meeting  the  control  objective  with  minimal 
amount  of  fuel  used  for  control. 


2.2.2  Two  dimensional  nonlinear  model 

Thermoacoustics  models  incorporating  effects  of  distributed  acoustics,  flame  dynamics  and  vorticity 
in  bluffbody  flameholder  annular  combustors  have  been  developed.  In  this  Section,  we  first  describe  a 
2D  nonlinear  model  without  vorticity  that  is  intended  for  liner  design.  Details  on  this  model  can  be 
found  in  [33],  where  first  a  3D  model  is  obtained  and  then  averaging  along  the  radial  direction  yields 
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a  2D  version  of  the  model.  The  linearized  potential  acoustics  for  this  model  are  described  by 
%  +  u  •  Vp  +  7  pST[Yf]S(x  -  xf)p  -  p0\v\U\2  •  V0 

,  .  -2a,  ,  dK\R*  _ 


-f  pqc A4>  + 


(7  - 1)?, 


^  +  u-V4>  = 


where  </>,p  denote  the  averaged  2D  acoustic  potential  and  pressure  respectively,  U  denotes  the  mean 
flow,  St  is  the  flame  speed,  a  function  of  mean  local  fuel  mass  fraction  Yjy  the  term  2  arises  on 

account  of  the  radial  boundary  conditions  and  model  the  effect  of  liner  (wall  normal  velocity),  and 
parameters  po>7,c  denote  fluid  density,  ratio  of  specific  heats,  speed  of  sound  respectively  and  the 
parameter  p  =  —  1).  q  represents  the  heat  release  perturbation,  which  arises  due  to  the  flame 

dynamics  modeled  by  G-equations  describing  the  motion  of  individual  flames 

dG 

+  [U  +  V<£]  •  VG  +  St\VG\  =  0,  (2.20) 

where  the  flame  speed  St  is  determined  from  the  solution  of  the  fuel  advection  equation,  written  in 
2D  as 

+  (U  +  V<^>)  •  VYj  =  0.  (2.21) 

The  boundary  conditions  arise  due  to  the  acoustic  boundary  condition  at  the  flameholder  walls 

V<£  ’  I  Flameholder  Walls  =  .  (2.22) 

and  the  inlet  fuel  profile 

Y?  = - — ^T-,  (2.23) 

/  /  _  1  V  \  (TT  ,  d(b\  ’  \ 


1  (P0+$)(U+$£Y 

where  p/Uf  denotes  the  fuel  mass  flux  and  U  denotes  the  axial  component  of  flow  velocity. 

In  [33],  we  present  the  above  model  for  describing  combustion  instabilities  in  annular  combustors. 
The  model  is  being  used  to  study  robust  linear  design,  where  the  control  input  acts  through  the 

boundary  liner  terms  ^  in  the  acoustic  model  above.  The  model  includes  a  model  for  mean  flow 

U  =  Up  +  Ue>  (2.24) 

where  the  potential  flows  Up  models  the  inflow  of  reactants  and  models  the  expansion  velocity 
created  because  of  burning,  arises  as  a  solution  to  the  continuity  equation 

1  Dp 

=  (2.25) 

together  with  appropriate  boundary  conditions  (see  [33]  for  details  on  the  model).  An  example  of  the 
above  model  with  only  a  single  flameholder  (see  Figure  2.3  for  schematic)  is  explicitly  constructed  and 
presented  in  the  report. 

Our  preliminary  attempt  at  extending  the  model  to  include  vortex  dynamics  is  to  use  the  decom¬ 
position  of  velocity  field  similar  to  equation  (2.24)  above  but  now  include  a  term  due  to  vorticity 

a^Hp  +  IU  +  ILn  (2.26) 

where  vortical  velocity  arises  as  a  solution  of  the  vorticity  equation 

1  1 

+  ([H  +  V0]  •  V)w  +  (V  •  U  -h  A(£)w  =  —Aw  +  ^2  Vp  X  VP.  (2.27) 

Rigorous  justification  of  the  above  extension  is  part  of  the  ongoing  research. 
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Figure  2.3:  Schematic  of  the  physical  problem. 


2.2.3  FLAVOR  -  A  vortex  model  for  reacting  flow  bluffbody  wake  dynamics 

In  our  paper  [36],  we  present  a  vortex  model  for  bluffbody  flameholder  stabilized  premixed  combustion. 
We  consider  the  physical  problem  of  non-premixed  combustion  stabilized  by  a  single  rectangular  bluff 
body  flameholder  of  height  h  in  a  channel  of  height  H  (see  Figure  2.3  for  a  schematic).  The  size  of 
both  the  bluff  body  and  the  channel  in  the  third  dimension  (z)  is  large  so  that  two-dimensionality 
is  assumed  to  apply.  The  Mach  number  is  low  and  both  reactants  (of  density  pu)  and  products  (of 
density  p &)  are  assumed  to  behave  as  ideal  gases.  The  combustion  time  scale  is  much  faster  than  that 
of  the  flow  and  the  reacting  field  is  assumed  to  be  approximated  by  two  flamesheets  anchored  at  the 
two  flameholder  lips. 

The  evolution  of  the  flow  is  governed  by  the  vorticity  and  continuity  equations 


dco 


(u  *  V)o;  +  (V  •  u)u> 


V  •  u 


±-Au,  +  ±VpxVp, 

_lDp 

p  Dt 


(2.28) 

(2.29) 


together  with  the  no-slip  and  impermiability  boundary  conditions  at  the  flameholder  walls  and  im- 
permiability  boundary  conditions  at  the  channel  walls.  A  G-equation  formulation  (see  [27])  is  used  to 
describe  the  flame  evolution  as 


_  +  ([M  +  ST(xf)n].V)G  =  0,  (2.30) 

where  the  flamesheet  is  described  by  the  (unidimensional)  connected  locus  of  the  points  xj  satisfying 
G(xj,t)  =  0  and  h  is  the  unit  normal  to  the  flame  oriented  into  reactants.  For  the  premixed  case 
considered  here,  the  flame  speed  is  modeled  to  be  its  stoichiometric  value,  while  retaining  the  effects 
of  curvature  (using  Markstein  length). 

The  Lagrangian  Vortex  Element  Method  (VEM)  in  the  form  that  accommodates  the  presence  of 
reaction  in  the  low  Mach  number  limit  [22]  is  used  to  reproduce  the  unsteady  flow.  The  Lagrangian 
implementation  of  the  flamesheet  evolution  is  performed  using  numerical  techniques  consistent  with 
the  VEM.  The  numerical  model  is  validated  against  non-reacting  experimental  data  for  shear  layers 
(results  included  in  [14])  simulated  by  considering  a  thin  bluff  body  with  a  velocity  difference  across 
it  -  and  traditional  bluff  body  flows  [44].  The  results  of  this  validation  is  presented  in  our  paper  [36]. 

Results  for  the  reacting  bluff  body  flow  indicate  a  shift  in  the  solution  from  the  Von  Karman  asym¬ 
metric  shedding  of  coherent  vortices  at  a  characteristic  frequency  witnessed  in  the  non-reacting  flow, 
to  a  rather  symmetric  shedding  that  is  not  dominated  by  any  single  frequency.  Figure  2.4  contrasts 
the  time-series  and  frequency  spectra,  based  on  the  v- velocity  signal  taken  on  the  centerline  half  bluff 
body  width  downstream  of  the  bluff  body  trailing  edge,  for  the  reacting  and  non-reacting  cases.  The 
reduction  in  unsteadiness  (with  respect  to  the  non-reacting  case)  seen  in  the  figure  is  consistent  with 
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Spectrum 


(a)  (b) 

Figure  2.4:  (a)  Time  series  and  (b)  spectral  plots  of  the  non-reacting  (blue)  and  reacting  (red)  v- 
velocity  signal  on  the  centerline  (y=0)  at  half  a  bluff  body  width  downstream  from  trailing  edge. 

the  experimentally  obtained  spectral  plots,  comparing  non-reacting  and  reacting  spectra,  presented  in 
[47].  Analysis  presented  in  our  paper  [36]  indicates  that  this  shift  is  mainly  due  to  the  dilatation  that 
accompanies  the  combustion  heat  release  while  baroclinic  vorticity  plays  a  supporting  but  secondary 
role.  The  dynamics  in  the  near  field  of  the  flameholder  (4-5  bluff  body  thickness  downstream)  is  dom¬ 
inated  by  the  vorticity  generated  at  the  bluff  body  walls.  The  dilatation  weakens  this  wall  generated 
vorticity  as  it  goes  through  the  flame  and  delays  the  entrainment  of  some  of  this  vorticity  into  the 
products  region  of  the  wake.  Both  effects  tend  to  diminish  the  interaction  of  the  opposite  signed  vor¬ 
ticity  emerging  from  the  boundary  layers  on  the  two  horizontal  bluff  body  walls,  thereby  diminishing 
the  possibility  of  a  Von  Karman  street.  Further  downstream  vorticity  generated  by  the  baroclinic 
torque  dominates  the  dynamics.  The  amount  of  the  downstream  baroclinic  vorticity  is  strongly  de¬ 
pendent  on  the  presence  of  wall  generated  vorticity.  The  latter  excites  the  flame  in  the  near  field 
thereby  enhancing  the  conditions  for  baroclinic  vorticity  generation.  Characteristic  simulation  results 
are  shown  in  Figure  2.5. 

2.2.4  Reduced  order  modeling  for  Control 

Fuel  control  is  a  viable  strategy  for  suppressing  the  combustion  instability  induced  pressure  oscillations. 
However,  the  application  of  fuel  control  to  bluffbody  environment  with  its  distributed  fluid  dynamics, 
combustion  and  acoustics  is  yet  to  be  done.  In  the  following  subsection,  we  present  the  past  research 
whose  explicit  aim  was  to  apply  reduced  order  modeling  for  the  purpose  of  investigating  control. 

A  primary  difficulty  in  applying  fuel  control  to  bluffbody  flameholders  arises  due  to  the  lack  of 
suitable  reduced  order  models  that  can  be  used  for  control  design.  In  some  of  the  past  research 
at  UTRC,  we  made  a  simplifying  assumption  of  neglecting  the  effect  of  vortical  dynamics  to  obtain 
reduced  order  models,  which  we  used  for  studying  fundamental  limitations  on  control  design.  The 
control  was  demonstrated  on  the  full  computational  model.  Below  we  provide  a  summary  of  this 
research;  see  [34]  for  details. 

We  consider  the  physical  problem  described  above  in  Section  2.2.3  and  in  figure  2.3.  The  primary 
difference  here  is  that  while  the  flow  is  simplified  (by  neglecting  vortex  dynamics),  combustion  now 
is  assumed  to  be  non-premixed  (ho  additional  equations  are  needed  to  describe  the  evolution  of  fuel). 
For  the  purpose  of  modeling  and  controlling  combustion  instability,  we  are  primarily  interested  in  the 
spatio-temporal  distribution  of  the  heat  release  response  q(x,t)  which  is  function  of  local  flame  speed 
and  the  flame  location 

q(x,t)  oc  St{xj)S{x  -£/).  (2.31) 

The  flame  dynamics  arise  due  to  (complicated)  interaction  between  the  fluid  dynamics  and  the 
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<b) 

Figure  2.5:  Flow  visualization  of  the  reacting  flow  with  h/H=0.16,  h/H=20,  Re=20000,  and  .  (a) 
Instantaneous  flame  location  (lines),  center  of  vortex  elements  (black/white  points  denoting  -/+  signed 
vorticity)  together  with  the  mean  streamwise  velocity  field  (shades),  (b).  Mean  flame  location  (lines) 
together  with  mean  vorticity  (shades)  -  the  color  axis  range  for  vorticity  plot  is  chopped  to  better  view 
the  weaker  bar o clinic  vorticity. 
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*  A  typical  flame  speed  function  of  fuel  mass  fraction 


Figure  2.6:  Plot  of  a  typical  flame  speed  function  St\Yj\  :  the  peak  corresponds  to  stoichiometric 
condition  where  the  flame  burns  with  maximal  flame  speed  So- 

inherent  kinematic  flame  motion.  The  2D  nature  of  the  problem,  boundary  conditions  and  bluffbody 
geometry  (where  there  are  no  translation  symmetries)  makes  reduced  order  modeling  for  the  purpose  of 
understanding  and  controlling  the  reacting  flow  problem  difficult.  We  make  a  simplifying  assumption 
of  ignoring  the  vortical  component  of  the  fluid  dynamics  in  the  problem,  but  retain  the  effects  due 
to  burning  and  fuel  actuation.  The  equations  of  flow  retain  the  continuity  equation  but  the  vorticity 
equation  is  now  replaced  with 

V  x  u  =  0,  (2.32) 

and  for  the  boundary  conditions,  only  the  impermiability  condition  is  retained  (the  no-slip  condition 
is  dropped).  The  G-equation  is  once  again  used  to  describe  the  flame  motion  but  now  for  the  non- 
premixed  case,  the  flame  speed 

Srixf)  =  Sr[Yf{xf)]  (2.33) 

is  a  function  of  the  local  fuel  mass  fraction  for  the  non-premixed  situation  considered  here.  Figure  2.6 
plots  a  typical  Sr[Yf]  as  a  function  of  Yj.  We  reserve  the  square  bracket  notation  St[-]  for  the  function 
to  distinguish  it  from  the  flame  speed  Sj\  The  local  fuel  mass  fraction  seen  at  the  flame  front  arises 
due  to  the  convection  of  the  fuel-air  mixture  in  the  duct 

T?  -  °-  <2-34> 

An  initial  upstream  profile 

Yf(x  =  0,y,t)=Yf0(y,t)  (2.35) 

provides  a  boundary  condition  at  the  entrance  of  each  of  the  half  channels  upstream  (on  either  side) 
of  the  bluff  body  (see  Figure  2.3).  This  initial  distribution  convects  with  the  flow  velocity  (because 
of  (2.34))  until  it  burns  at  the  flame  front.  Here,  we  assume  that  the  fuel  distribution  is  such  that  a 
lean  condition  (where  the  fuel  mass  fraction  Yp  is  less  than  the  stoichiometric  value-  see  Figure  2.6) 
always  applies.  As  a  result,  all  of  the  fuel  is  burnt  at  the  flame  and  the  flame  provides  an  appropriate 
second  boundary  condition  for  (2.34). 
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Figure  2.7:  Plot  showing  shallowing  of  flame  angle  as  density  ratio  increases:  the  difference  between 
numerical  model  (~)  and  reduced  order  model  (dashed)  is  less  than  0.5  degrees  even  for  large  density 
ratios. 

We  derive  a  reduced  order  model  for  the  individual  flame  motion  as 

!  +  {Sr[l?(*,*- £)]+!#<  +  (2.36) 

€(v  =  0,t)  =  0,  (2.37) 

£(»,*  =  0)  =  £°,  (2.38) 

a  hyperbolic  initial  boundary  value  problem  (IB VP)  -  here  (£(y,t),?/)  denote  the  co-ordinate  of  the 
flaune  location.  Here  /i  =  (^  —  1),  where  pu  is  the  density  of  the  unburnt  reactants  and  p^  is  the 
density  of  the  products. 

We  also  validated  the  reduced  order  model  against  the  simulation  results  of  a  CFD  numerical 
model.  The  reduced  order  model  was  shown  to  accurately  reproduce  both  the  flame  shapes  and  flame 
angles  as  a  function  of  density  ratio  parameter  p  for  premixed  flames  (see  Figure  2.7)  and  accurately 
predict  the  flame  blowoff  seen  in  the  computational  model  for  the  non-premixed  cause  as  the  fuel 
penetration  into  the  cross-stream  increases  (see  Figure  2.8). 

We  also  used  the  reduced  order  model  to  define  and  study  a  (fuel  actuation  based)  control  problem 
of  tracking  the  heat  released  due  to  flame  motion  against  a  prescribed  reference  signal  (for  the  full 
simulation).  The  problem  is  motivated  by  the  problem  of  fuel  control  of  combustion  instability.  In 
order  to  completely  address  the  issue  of  control  authority  needed  to  quench  pressure  oscillations  in 
combustor  to  an  acceptable  level,  one  needs  to  investigate  the  coupled  acoustic-flame  system  model 
[32].  In  the  paper,  however,  we  concentrated  on  controlling  the  distributed  heat  release  response  by 
manipulating  the  inlet  fuel  profile  (see  Figure  2.9  for  a  schematic  of  the  control  problem).  A  heuristic 
interpretation  for  considering  this  piece  of  the  problem  is  provided  in  [34].  For  the  control  problem 
described  in  the  figure,  we  proposed  an  optimal  control  strategy  which  was  then  implemented  on  the 
computational  model  to  track  the  heat  released  against  a  prescribed  reference  signal.  Figure  2.10 
provides  a  comparison  of  the  tracking  signal  against  the  heat  release  obtained  from  the  simulation. 
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Figure  2.8:  Flame  angle  as  fuel  center  ymi(i  (defines  the  inlet  fuel  profile  as  Yj(y)  =  JS-  e  <3<d)  ) 
varies:  model  captures  the  critical  value  of  flame  blow-off  as  well  as  flame  angles  for  robust  flame 
before  blow-off. 


(control  input) 


Figure  2.9:  Schematic  of  the  control  problem. 
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Desired(-)  and  numerical  model(— ) 


Figure  2.10:  Results  of  the  implementing  the  tracking  control  on  the  numerical  model. 

Analysis  of  the  control  was  provided  in  the  paper  and  suggested  in  a  fundamental  limitation  associated 
with  ON-OFF  control  of  heat  release.  This  limitation  is  suggested  in  Figure  2.10  as  discrepancy  be¬ 
tween  the  desired  and  the  controlled  heat  release  response.  The  limitation  arises  due  to  the  hyperbolic 
nature  of  the  problem  which  leads  to  solution  discontinuities  with  ON-OFF  control;  complete  analysis 
using  method  of  characteristics  is  provided  in  [34]. 

2.2.5  Reduced  order  modeling  for  Dynamic  range  analysis 

The  occurrence  of  combustion  instability  in  augmentors  using  bluffbody  flameholders  is  a  function 
of  the  operating  conditions  [8],  Dynamic  range  analysis  is  important  not  only  to  map  the  stability 
boundaries  (where  instabilities  occur)  and  predict  post  instability  response  (for  instance,  the  pressure 
amplitudes)  but  also  to  identify  the  important  physical  mechanisms  that  lead  to  the  instabilities  as  a 
function  of  parameters  in  the  problem.  There  are  two  challenges  associated  with  meeting  the  objective 

1.  The  computational  cost  associated  with  mapping  out  all  the  different  regimes  of  physical  phe¬ 
nomenon  is  prohibitive, 

2.  Even  if  a  transition  of  physical  phenomenon  is  seen  in  the  computational  model,  an  understand¬ 
ing  of  key  physical  mechanisms  underlying  this  transition  is  difficult  to  obtain  because  of  the 
distributed  and  nonlinear  nature  of  the  dynamics. 

For  the  purpose  of  reducing  the  computational  burden,  we  are  interested  in  model  reduction  approaches 
that  can  model  the  computationally  expensive  elements  (such  as  vortex  dynamics)  in  reduced  order 
fashion.  In  the  following  subsection,  we  summarize  one  such  approach  -  a  Galerkin  based  reduced 
order  modeling  framework  -  used  to  understand  the  suppression  of  Von  Karman  vortex  shedding  for 
compressible  reacting  flows.  This  ongoing  work  is  as  yet  unpublished,  but  available  as  a  preprint  [35]. 

The  wake  dynamics  of  a  bluffbody  stabilized  reacting  flows  are  very  different  from  their  non¬ 
reacting  counterpart.  In  particular,  while  a  Von  Karman  asymmetric  shedding  of  coherent  vortices 
at  a  characteristic  frequency  is  witnessed  in  the  non-reacting  flow,  reacting  flows  exhibit  a  rather 
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symmetric  shedding  that  is  not  dominated  by  any  single  frequency  (see  our  recent  computational 
paper  [36],  also  the  review  paper  [10],  computational  papers  of  [9]  and  [41]  and  experiments  of  [47]). 
A  variety  of  explanations  have  been  proposed  in  the  literature,  the  majority  of  them  pointing  to 
the  vorticity  generated  by  the  flame  via  the  baro clinic  mechanism  as  the  primary  reason  behind 
the  exhibited  shift  in  flow  behavior.  For  example,  [10]  attributes  the  shift  to  the  combination  of  two 
possible  mechanisms:  (i)  the  dampening  of  the  vorticity  due  to  the  increased  kinematic  viscosity  of  the 
reacted  fluid  and,  (ii)  the  generation  of  baroclinic  vorticity  that  is  of  opposite  sign  to  the  flameholder 
generated  vorticity  and  tends  to  nullify  the  effects  of  the  latter.  Menon  and  co-  workers  also  point 
to  the  baroclinic  vorticity  generation  as  the  main  mechanism  that  leads  to  the  shift  [9].  In  our  own 
paper  [36],  we  argue  the  importance  of  exothermic  effects  in  suppressing  the  Von  Karman  shedding 
observed  in  cold  bluffbody  flow. 

In  our  more  recent  work,  we  use  Galerkin  based  reduced  order  models  for  investigating  the  effect 
of  exothermicity  on  reacting  bluffbody  flows.  In  particular,  we  are  interested  in  explaining  (within 
the  reduced  order  models)  the  suppression  of  Von  Karman  shedding  in  the  presence  of  burning.  In 
order  to  obtain  model  reduction,  our  approach  is  to  reduce  the  order  of  the  fluid  dynamics  using 
(Noack-Tadmor  [21])  inspired)  POD  based  Galerkin  model.  To  obtain  the  Galerkin  model,  we  begin 
with  the  compressible  form  of  vorticity  equation 

+  (w  *  V)o/  +  (V  •  u)co  =  —  Acj  +  Vp  x  Vp.  (2.39) 

We  make  a  simplifying  assumption  of  ignoring  the  effect  of  baroclinicity, 


-Vp  x  Vp  —  0, 
P 


(2.40) 


which  is  motivated  by  some  of  the  results  summarized  in  [36],  where  we  have  argued  the  importance 
of  exothermicity  (V  •  u  >  0)  in  obtaining  suppression. 

In  [35],  we  show  that  the  reduced  order  Galerkin  model  (with  three  POD  modes)  has  the  following 
structure 
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where  (waj^i?^)  are  the  modal  coefficients,  the  first  term  on  the  right  hand  side  of  equation  (2.41) 
models  the  non-reacting  wake  flow  dynamics  (see  [21]  for  their  physical  significance)  and  the  second 
and  third  terms  model  the  effect  of  burning  (recall  p,  =  (^  —  1)).  The  second  term  models  the  effect  of 
exothermicity  captured  by  the  compressible  term  (V  •  u)u>  in  the  vorticity  equation  (2.39).  The  third 
term  models  the  effect  of  convection  because  of  the  dilatation  flow  (arises  as  a  component  of  the  term 
(u  ■  V)o>  in  the  vorticity  equation  (2.39))  which  causes  the  vorticity  to  convect  differently  than  for  a 
non-reacting  flow. 

In  the  absence  of  burning,  the  3  mode  Galerkin  model  reproduces  the  sequence  of  bifurcations 
from  fully  attached  steady  flow  to  a  steady  flow  with  symmetric  recirculation  regions  in  the  wake 
to  an  assymetric  Von  Karman  shedding  solution  (see  Figure  2.11).  In  the  reduced  order  model,  the 
effect  of  the  burning  is  to  laminarize  the  flow  by  moving  the  Hopf  bifurcation  point  corresponding  to 
the  appearance  of  Von  Karman  shedding  solution  to  higher  Reynolds  numbers.  Furthermore,  beyond 
a  critical  value  of  parameter  (p  «  1.4),  the  shedding  solution  dissapears  for  all  values  of  Reynolds 
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Figure  2.11:  Bifurcation  diagram  for  the  Galerkin  model 


numbers  (see  Figure  2.12  for  a  locus  of  Hopf  bifurcation  point  as  a  function  of  the  two  parameters  - 
Reynolds  number  Re  and  /i).  The  value  of  /i  =  1.4  is  in  a  nice  agreement  with  our  reported  value  for 
disappearance  of  Von  Karman  shedding  solution  for  high  Reynolds  number  flows  [36]. 

2.3  Passive  control  of  thermoacoustic  instabilities  by  symmetry¬ 
breaking 

In  paper  by  Hagen  and  Banaszuk,  CDC  2004,  we  presented  a  thermo-acoustic  model  on  a  cylindrical, 
or  annular,  geometry,  capable  of  modeling  instabilities  of  tangential  acoustic  modes.  The  model 
accounts  for  non-uniform  density,  damping,  rotational  flow,  and  heat-release  coupling.  It  is  shown  that 
deliberately  introducing  spatial  variations  in  some  quantities  has  a  similar  effect  to  adding  damping 
to  the  system.  The  effects  of  these  symmetry-breaking  conceptes  are  evaluated  on  the  model  through 
linear  analysis  and  the  net  amount  of  additional  damping  is  computed.  We  showed  how  various 
symmetry-breaking  concepts  are  robust  with  respect  to  the  uncertainty  in  the  model  parameters  and 
we  examined  propagation  of  uncertainty  with  respect  to  a  measure  of  uncertainty  recently  defined  by 
I.  Mezic. 

2.4  Background  noise  effect  on  combustor  stability 

Paper  by  Lieuwen  and  Banaszuk,  Journal  of  Propulsion  and  Power  2004,  considers  the  effects  of  back¬ 
ground  turbulent  fluctuations  upon  a  combustor’s  stability  boundaries.  Inherent  turbulent  fluctuations 
act  as  both  additive  and  parametric  excitation  sources  to  acoustic  waves  in  combustors.  While  ad¬ 
ditive  noise  sources  exert  primarily  quantitative  effects  upon  combustor  oscillations,  parametric  noise 
sources  can  exert  qualitative  impacts  upon  its  dynamics;  particularly  of  interest  here  is  their  ability 
to  destabilize  a  system  that  is  stable  in  the  absence  of  these  noise  sources.  The  significance  of  these 
parametric  noise  sources  increases  with  increased  background  noise  levels  and,  thus,  may  play  more  of 
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Figure  2.12:  Locus  of  Hopf  bifurcation  point  in  Re-/*  plane:  The  Hopf  point  denoting  the  Reynolds 
number  where  the  Von  Karman  vortex  shedding  solution  begins  is  moved  out  to  its  inviscid  limit  as 
/x  increases  to  «  1.4. 

a  role  in  realistic,  high  Reynolds  number  systems  than  experiments  on  simplified,  lab  scale  combustors 
might  suggest.  The  objective  of  this  paper  is  to  determine  whether  and/or  when  these  effects  might 
be  significant.  The  analysis  considers  the  effects  of  fluctuations  in  damping  rate,  frequency  and  com¬ 
bustion  response.  It  is  found  that  the  effects  of  noisy  damping  and  frequency  upon  the  combustor’s 
stability  limits  is  relatively  small,  at  least  for  the  fluctuation  intensities  estimated  here.  The  effects 
of  a  noisy  combustion  response,  particularly  of  a  fluctuating  time  delay  between  flow  and  heat  release 
perturbations,  can  be  quite  significant,  however,  in  some  cases  for  turbulence  intensities  as  low  as  5- 
10%.  These  results  suggest  that  deterministic  stability  models  calibrated  on  low  turbulence  intensity, 
lab  scale  combustors  may  not  adequately  describe  the  stability  limits  of  realistic,  highly  turbulent 
combustors. 


2.5  Modeling,  Analysis,  and  Control  of  Flutter  in  Turbomachinery 

A  linear  framework  for  control  of  wave  phenomena  on  annular  domain  was  established.  The  flutter 
control  problem  was  used  to  motivate  the  study.  However,  the  framework  applies  to  control  of  general 
wave  phenomena  on  annular  domain,  such  as  rotating  stall  and  thermoacoustic  instabilities. 

In  papers  by  Banaszuk  et  a/.,  IFAC  2002,  CDC  2002,  AIAA  2003,  we  described  a  method  for 
controlling  fan  or  compressor  blades  flutter  in  gas  turbine  engines  and  its  experimental  demonstration. 
The  experimental  implementation  of  active  flutter  control  on  a  sub-scale  fan  rig  consisted  of  an  array 
of  audio  speaker-powered  volumetric  sources  connected  to  the  flow  path  equally  spaced  along  the 
circumference  axially  located  between  an  experimental  fan  and  its  exit  guide  vanes  (stator).  Blade 
arrival  time  detectors  based  on  eddy  current  sensors  placed  at  the  leading  end  of  the  blade-tip  line 
were  used  to  generate  real-time  blade  deflection  signals.  An  observer  was  used  to  reconstruct  the 
flutter  modes.  A  pole-placement  controller  was  used  to  generate  the  speaker  command  signals.  The 
control  system  was  able  to  add  significant  amount  of  damping  to  three  modes  of  flutter.  Damping 
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augmentation  was  an  order  of  magnitude  larger  than  the  intrinsic  aeromechanical  damping  of  the 
modes  on  the  operating  line  of  the  fan. 

Blade  failures  due  to  flow  induced  vibrations  axe  a  long  standing,  endemic  problem  for  the  turbo¬ 
machinery  industry.  Flutter  and  resonant  stress  fundamentally  constrain  the  design  and  operation  of 
gas  turbine  engines.  Ensuring  aeromechanical  operability  often  requires  compromises  in  turbomachine 
efficiency,  performance  and  cost  and  can  result  in  development  delays  and  increased  mainatinance 
costs. 

We  describe  a  method  for  controlling  fan  and  compressor  balde  flutter  in  gas  turbine  engines,  as 
well  as  a  particular  implementation  of  this  approach  demonstrated  in  experiments  in  a  transsonic  fan 
rig  operating  at  9000  RPM.  The  limitations  in  operability  of  a  turbomachine  due  to  flutter  can  be 
overcome  by  adding  damping  to  the  dominant  aeromechanical  modes. 

We  model  the  dynamics  of  blade  rows  in  turbo-machinery  as  similar  to  those  of  a  flexible  disk. 
Aeromechanical  modes  form  travelling  waves  as  seen  by  the  rotor.  This  means  that  when  viewed  from 
the  rotating  frame  the  peak  of  deflection  appears  to  travel  around  the  disk.  The  deflection  of  the  disk 
at  a  given  point  on  the  fixed  frame  along  the  circumference  of  the  blade-row  can  be  decomposed  into 
sinusoids  of  frequencies  separated  by  integer  multiples  of  the  rotor  frequency.  At  any  fixed  point  in 
time,  the  deflection  of  the  disk  can  also  be  decomposed  into  sine-waves  function  of  the  angular  position 
around  the  rotor.  Therefore  each  aeromechanical  mode  has  a  characteristic  shape  and  a  characteristic 
frequency.  Each  of  these  modes  can  lose  stability  as  operating  conditions  change.  The  objective  of  a 
flutter  control  is  to  enhance  the  region  of  stable  operation  by  adding  damping  to  the  aeromechanical 
modes. 

The  preferred  sensing  scheme  uses  a  proximity  sensor  on  the  casing  to  determine  a  time  blade 
arrival  time.  From  the  blade  arrival  time  one  can  estimate  the  blade  deflection.  The  blade  deflection 
will  reflect  the  superposition  of  all  the  aeromechanical  modes.  However,  due  to  their  separation  in 
frequency,  the  modal  content  can  easily  be  decomposed. 

The  actuation  approach  is  to  place  volumetric  sources  aft  of  the  blade-row  to  modulate  the  back 
pressure  and  mass  flow  as  a  function  of  angular  position  and  time  resulting  in  unsteady  loading  of  the 
blades.  This  in  turn  modifies  the  blade  lift,  generating  the  desired  commanded  force  on  the  blades.  By 
arranging  an  array  of  such  actuators  around  the  circumference  one  can  create  a  pattern  of  forces  on 
the  blades.  These  patterns  can  form  traveling  waves  that  have  the  spatial  shape  of  the  aeromechanical 
modes.  The  experimental  implementation  of  active  flutter  control  on  a  fan  rig  presented  in  this  paper 
consisted  of  an  array  of  audio  speaker-powered  volumetric  sources  connected  to  the  flow  path  equally 
spaced  along  the  circumference  axially  located  between  the  fan  rotor  blades  and  its  exit  stator  guide 
vanes. 

In  experiments  on  a  fan  rig  a  linear  observer  was  implemented  to  estimate  the  aeromechanical 
modal  content  of  the  blade  row.  This  approach  required  a  linear  model  for  the  dynamics  of  interest. 
Such  a  model  was  obtained  for  each  aeromechanical  mode  by  running  swept-sine  experiments  in  the 
system  and  measuring  the  complex  ratio  between  the  modal  forcing  function  and  the  blade  deflection 
at  a  point  in  the  fixed  frame.  Then  a  low  order  (typically  second  order)  state-space  linear  system 
was  fit  to  the  frequency  response  data.  The  observer  was  designed  based  on  the  aggregate  of  all  these 
state-space  blocks.  An  observer-based  pole-placing  technique  was  used  to  design  a  linear  control  law 
to  add  the  desired  amount  of  damping  to  flutter  modes.  The  control  system  was  able  to  add  damping 
to  three  flutter  modes.  The  damping  augmentation  achieved  was  an  order  of  magnitude  larger  than 
the  intrinsic  aeromechanical  damping  of  the  flutter  modes  on  the  most  stable  operating  condition  of 
the  fan. 

Experiments  described  in  this  section  were  conducted  under  internal  funding.  Modeling  and  anal¬ 
ysis  were  partially  funded  under  previous  AFOSR  grants.  The  publications  Banaszuk  et  a/.,  IFAC 
2002,  CDC  2002,  AIAA  2003  of  the  results  were  funded  under  current  AFOSR  grant. 
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2.5.1  Notation 


N  -  number  of  blades. 

n  -  index  of  a  flutter  mode,  n  =  . . . ,  —2,  —1, 0, 1, 2, _ 

— Cn  -  real  part  of  the  n-th  flutter  mode  pole. 

unr  -  imaginary  part  of  the  n-th  flutter  mode  pole,  circular  (pseudo)  frequency  of  the  n-th  flutter 
mode  in  the  rotating  frame. 

tons  -  imaginary  part  of  the  n-th  flutter  mode  pole,  circular  frequency  of  the  n-th  flutter  mode  in  the 
stationary  frame. 

£n  -  damping  coefficient  of  the  n-th  flutter  mode,  £n  := 

Sn  -  Logarithmic  decrement  of  the  n-th  flutter  mode,  8n  :=  M  =. 

9r  -  angle  in  the  rotating  frame. 

9S  -  angle  in  the  stationary  frame. 

anr(t,9r)  -  blade  deflection  angle  at  time  t  at  angle  9r  (in  the  rotating  frame). 

&ns{t, 9S)  -  blade  deflection  angle  at  time  t  at  angle  9S  (in  the  stationary  frame). 

(vr  -  circular  rotor  frequency. 

8so  -  angle  between  the  fixed  reference  points  on  the  rotor  and  the  stator  at  time  t  =  0. 

(*)n  -  n-th  spacial  Fourier  coefficient. 

(:)  -  temporal  Fourier  transform. 


2.5.2  Flutter  models 


For  an  integer  n  (positive,  zero,  or  negative)  we  model  the  n-th  flutter  mode,  or  n-th  nodal  diameter 
flutter  mode,  as  a  travelling  wave  in  which  all  blades  are  oscillating  harmonically  with  a  constant 
phase  angle  0n  :=  relative  to  each  other  [18]. 

Let  0r  denote  the  angle  measured  relative  to  a  fixed  point  on  the  rotor  in  the  direction  of  the 
rotation.  Assume  that  we  have  continuum  of  blades  and  there  is  no  external  forcing.  We  postulate 
that  the  n-th  nodal  component  of  the  blade  deflection  at  angle  9r  at  time  t  is  given  by  the  formula 

^nr(^)^r)  ”  COS (u)nrt  —  Tl6r  +  (f>nr )  (2.42) 

where  —  £n  and  unr  are,  respectively,  the  real  and  imaginary  part  of  the  n-th  flutter  mode  pole.  Note 
that  Ljnr  is  also  the  (pseudo)  frequency  of  the  n-th  flutter  mode  in  the  rotating  frame,  and  An  and  <f>nr 
are  the  initial  magnitude  and  phase  angle  of  the  n-th  flutter  mode.  The  damping  of  the  n-th  flutter 
mode  is  usually  described  by  one  of  two  coefficients:  the  damping  coefficient  £n  :=  r  &  or  the 
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logarithmic  decrement  5n  :=  2^-^-  =  27 r  .  Note  that: 

(1)  The  m-th  blade  is  moving  according  to  equation  (2.42)  with  the  corresponding  angle  6r  =  + 

where  0i  is  the  position  of  the  first  blade  relative  to  the  fixed  reference  point  on  the  rotor,  m  = 
1,2,... ,AT. 

(2)  For  a  fixed  time  t  and  n/0  the  blade  deflection  anr{t,6r)  considered  as  a  function  of  the  angle 
0T  has  a  sinusoidal  shape  with  |n|  nodes.  For  n  =  0  and  a  fixed  time  t  the  deflection  is  the  same  for 
each  blade. 

(3)  For  =  0  and  n  ^  0  the  blade  deflection  aw(£>0r)  is  a  wave  with  a  fixed  sinusoidal  shape 
travelling  around  the  annulus.  The  speed  and  the  direction  of  rotation  can  be  obtained  by  considering 
movement  in  time  of  the  angle  corresponding  to  one  of  the  peaks  of  the  wave.  For  instance,  the  first 
peak  is  obtained  by  solving  the  equation  wnrt  —  n0r  4-  (j>nr  =  f  for  6r.  We  have 


9r  —  {u)nrt  “h  <f>n 
n 
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Therefore,  the  speed  of  the  wave  is  ^  and  the  direction  is  positive  (the  same  as  the  direction  of 
rotation  of  the  rotor)  for  n  >  0  and  negative  (the  opposite  to  the  rotor’s  rotation  direction)  for  n  <  0. 
We  call  the  flutter  modes  travelling  in  the  same  direction  as  the  rotor  the  forward  travelling  modes 
and  the  ones  travelling  in  the  direction  opposite  to  the  rotor’s  direction  the  backward  travelling  modes. 
(4)  For  a  fixed  angle  9r ,  the  blade  deflection  anr(t,9r)  considered  as  a  function  of  time  represents  a 
response  of  a  damped  oscillator,  i.e.,  a  second  order  system  with  poles  —  £n  +  iunr  and  —  £n  —  iunr. 
Note  that  each  particular  blade  oscillates  with  frequency  a>nr,  which  is  n  times  bigger  than  the  fre¬ 
quency  of  the  corresponding  travelling  wave. 

Now  we  express  the  motion  of  a  blade  due  to  a  particular  flutter  mode  as  measured  at  an  arbitrary 
angle  on  the  stator. 

Let  u>r  denote  the  circular  rotor  frequency.  Fix  a  reference  point  on  the  stator.  The  angles  in  the 
stationary  frame  will  be  measured  relative  to  this  point  with  positive  direction  corresponding  to  the 
rotor’s  rotation  direction.  Let  6S o  denote  the  angle  at  which  the  reference  point  on  the  stator  is  seen 
from  the  reference  point  on  the  rotor  at  time  t  =  0.  Then,  for  an  arbitrary  time  t,  a  fixed  angle  0S  on 
the  stator  is  related  to  the  corresponding  point  on  the  rotor  0r  (measured  in  the  rotating  frame)  by 
the  formula  6r  =  6S  -b  6S o  —  u>rt.  Therefore,  the  deflection  of  the  blade  passing  a  fixed  angle  0S  on  the 
stator  at  time  t  is  given  by  the  formula 

9$)  =  &nr(f'')9s  +  9S o  COrt )  =  . 

Ane~ Cnt  cos((o;nr  +  mor)t  -  n9s  +  <f>ns)  '  *  ' 

where  (f>ns  :=  <pnr  —  n9s o  is  the  initial  phase  of  the  n-th  mode  in  the  stationary  frame. 

Note  that  for  =  0  and  n  ^  0  the  blade  deflection  ans(£,  9S)  in  the  stationary  frame  is  a  wave  with 
a  fixed  sinusoidal  shape  travelling  around  the  rotor.  In  particular,  a  single  blade  vibration  frequency 
in  the  stationary  frame  is  uns  :=  ujnr  +  no y. 

The  velocity  of  the  rotation  of  the  wave  can  be  obtained  in  a  similar  manner  as  in  the  rotating 
frame  case.  In  particular,  the  velocity  of  the  wave  corresponding  n-th  flutter  mode  in  the  stationary 
frame  is  ur  +  Let  us  recall  that  the  latter  is  the  velocity  at  which  a  fixed  point  on  the  graph 
of  the  blade  deflection  as  a  function  of  angle  (say,  a  peak)  is  travelling  around  the  annulus  at  the 
stationary  frame.  This  velocity  should  not  be  confused  with  an  individual  blade  velocity  due  to  n-th 
flutter  modes,  i.e.,  n;ns,  which  is  n  times  bigger. 

In  the  sequel  we  are  going  to  use  the  stationary  frame  only.  Therefore,  we  will  often  skip  the 
subscript  s  and  use  9  to  denote  the  angles  measured  in  the  stationary  frame. 

Since  at  a  fixed  time  the  flutter  modes  and  the  corresponding  forcing  functions  have  a  fixed 
sinusoidal  shape,  they  can  be  represented  via  their  spatial  Fourier  coefficients  (SFCs).  One  complex 
Fourier  coefficient  can  be  used  to  describe  a  single  sinusoidal  travelling  wave.  A  general  travelling  wave 
with  n-th  nodal  spatial  shape  and  a  temporal  frequency  ljq  has  the  form  /n(t,  9)  :=  Fn(t)  cos(u;o t~n9+ 
4>)  =  Fn(t)  cos(u;o£  +  <£)  cos(n0)  +Fn(t)  sin(u;ot  +  <£)  sin(n0).  The  corresponding  SFC  is,  for  n  ^  0,  equal 
to  /„(<)  :=  fn{t,0)ejn6de.  One  has  /»(*)  =  if  = 

^Fn(t)  /027r(eJ(w°t+^)  +  e-}(uot-2n8+4,)-)d0 

=  Fn(t){cos(w0t  +  <j>)  +  j  sin {w0t  +  <£)).  For  n  =  0,  one  has  f0(t)  :=  ±  /02ir  f0(t,  6)d9.  Thus,  f0(t)  = 
fo2*  Fo(t)  cos(u>ot  +  d>)d,0  =  Fo(t)  cos ( wot  +  <t>)  =  fo(t,  6),  for  all  6. 

To  reconstruct  a  wave  from  its  SFC  one  can  use  the  inverse  spatial  Fourier  tranform 

/»(*, 9)  =  Re(fn(t)*e?n6)  -  Re(fn(t)e^n9),  (2.45) 

where  (■)*  stands  for  the  complex  conjugation. 

Observe  that  for  n  0: 

(1)  The  magnitude  and  phase  of  the  complex  number  representing  the  spatial  Fourier  coefficient 
of  the  wave  /n(t,  0)  are  the  same  as  magnitude  and  phase  of  the  wave. 
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(2)  The  real  and  imaginary  part  of  the  spatial  Fourier  coefficient  of  the  wave  /n(£,  0)  are  the  Fourier 
series  coefficients  of  /n(£,0),  i.e.,  the  coefficients  of  /n(£,  0)  represented  as  a  linear  combination  of 
cos(n0)  and  sin(n0),  respectively. 

Assume  that  the  magnitude  and  phase  of  the  wave  /„(£,  0)  are  constant  in  time  with  Fn(t)  :=  Fn,  for 
some  n/0.  Then  /n(t,  0),  and  hence  /„(£),  is  a  periodic  function  of  t  and  one  can  define  the  temporal 

Fourier  transform  of  the  spatial  Fourier  coefficient  of  the  wave  /n(t,  0)  fn(jtv)  ;=  fn(t)e~^u,tdt  := 
f^Fn(t)ertWot+&e~ju)tdt  :=  Fn(t)e^6( u  —  uo),  where  £(•)  stands  for  the  delta  operator.  Thus,  the 
travelling  waves  with  the  temporal  frequency  u>o  can  be  recognized  in  the  (temporal)  frequency  domain 
as  “spikes”  at  one  single  frequency  uq.  Spikes  at  positive  frequencies  represent  the  forward  travelling 
waves,  whereas  the  spikes  at  negative  frequencies  represent  the  backward  travelling  waves. 

The  case  n  =  0  is  different.  As  we  have  noticed  before,  the  spatial  Fourier  coefficient  /o(£) 
of  the  function  /o(£,  0)  coincides  with  the  function  /o(£,  0)  itself.  Its  temporal  Fourier  transform  is 

/oCH  ==  J-oc,  fo(t)e~jutdt  =  f^FolQlieX^+V  +  e~^  t+^)e~iutdt 

=  a»o) +e~^8(u)+u>o)).  One  observes  that  the  temporal  Fourier  transform  of  the  spatial 

Fourier  coefficient  of  the  function  /o(£,  0)  has  two  “spikes”:  one  at  u>o  and  the  other  at  — wo- 

While  the  flutter  modes  for  n  ^  0  are  represented  by  travelling  waves,  they  can  be  excited  by 
forcing  inputs  that  are  either  travelling  waves  of  the  form  /n(t,0)  :=  Fn  cos(a;o£  —  nd  +  <f>)  or  by  the 
standing  waves  of  the  form 

/n(*>  0)  :=  Fn  cos(u)ot  +  <f>)  cos (n0).  (2.46) 

This  is  due  to  the  fact  that  a  standing  wave  can  be  represented  as  linear  combination  of  two  travelling 
waves:  Fn  cos(cj0£  +  <f>)  cos (n0)  =  ^i?Tn(cos(a;o£  —  nO  ■+■  4>)  +  cos(cl>0£  +  nd  +  <£)). 

The  temporal  Fourier  transform  of  the  standing  wave  (2.46)  is  fn(jui)  :=  = 

roo^n(^(ei(WnH^4e'J^^))e^ 

=  —urn)-he  ^6(co-hcjn)).  Note  that  the  latter  formula  is  valid  for  all  integers  n,  including 

n  =  0. 


2.5.3  Flutter  models  with  control 

We  assume  that  we  have  continuum  of  actuators  around  the  stator  that  influence  flutter  modes.  We 
will  control  the  n-th  flutter  mode  with  a  control  function  u(t,  0)  that,  as  a  function  of  angle,  has 
the  same  shape  as  the  n-th  flutter  mode  wave.  The  control  magnitude  and  phase  will  be  chosen 
appropriately  as  functions  of  the  measured  (or  reconstructed  using  an  observer)  magnitude  and  phase 
of  the  n-th  flutter  mode,  the  angle  in  the  stationary  frame,  previously  denoted  by  6S.)  Similarly,  for 
the  identification  purposes,  one  can  force  the  n-th  flutter  mode  with  a  using  a  wave  of  the  with  a 
constant  magnitude  and  phase.  More  precisely,  assume  that  the  control  input  forcing  function  for  the 
n-th  mode  is  a  travelling  wave  having  some  temporal  frequency  ujq  and  having  the  same  shape  as  the 
n-th  flutter  mode: 

un(t ,  0)  =  Un  cos(cJot  +  <f>nu  —  nd)  = 

Un  cos(uot  +  (f>nu)  cos(n0)+  (2.47) 

Un  sin (ojQt  4-  <j)nu)  sin(n0), 

for  some  constant  Un  and  (f)nu .  The  SFC  of  this  forcing  function  is  un(t)  =  Une^UQt^nu\  The 
corresponding  temporal  Fourier  transform  is  un{jw)  =  Une^^nu 8(u  —  c^o). 

We  also  assume  that  the  steady-state  n-th  flutter  mode  component  of  the  blade  deflection  response 
to  the  n-th  nodal  forcing  of  the  form  (2.47)  is  a  travelling  wave  with  the  same  spatial  shape  and 
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temporal  frequency,  possibly  shifted  in  phase  by  some  angle  cf)n  relative  to  the  forcing  function: 

an(£,  0)  =  An  cos(a;o*  -  nO  +  <£n )  = 

An  cos(a>ot  +  0n)  cos(n0)+  (2.48) 

An  sin(cJo^  4*  <f>n)  sin(n0) 

for  some  constant  An  and  cj)n  that,  for  fixed  Un  and  </>nu,  are  functions  of  u)q. 

The  SFC  of  the  ra-th  component  of  the  blade  deflection  is  an{t)  =  The  temporal 

Fourier  transform  of  the  SFC  of  the  n-th  component  of  the  blade  deflection  is  an(ju)  =  Ane^nS(u  — 
wo).  We  assume  that  we  measure  the  blade  displacement  at  finite  number  of  locations  on  the  stator. 
(This  is  going  to  be  accomplished  with  eddy  current  sensors.)  At  a  fixed  angle  9y  the  measured  blade 
dispalcement  is  going  to  be 

VnOyit)  :=  Otn(t^0y)  = 

An  cos(wot  —  nOy  +  (f>n)  ,  , 

~  An  cos (v0t  +  <f>n)  cos(n0y)+  '  ' 

An  sin(w0t  -f  (j>n)  sin(n0y). 

The  temporal  Fourier  transform  of  the  output  function  is  yn(jw)  =  ^-(e^^n~n0y^S(uj— wo)+e“J^n_n0»)j(w4- 

«*)). 

Now  we  present  dynamic  system  models  for  the  evolution  of  the  n-th  flutter  mode  subject  to 
control.  The  description  adapts  an  approach  to  model  rotating  stall  from  [42], 

One  can  obtain  a  low  order  model  describing  the  dynamics  of  the  n-th  flutter  suitable  for  control 
purposes  in  the  following  three  steps. 

1.  Conduct  an  experiment  to  obtain  the  transfer  function  between  the  n-th  SFC  of  the  forcing 
function  given  by  (2.47)  and  the  corresponding  n-th  SFC  of  the  blade  deflection  function  given  by 
(2.48). 

2.  Fit  a  low-order  transfer  function  to  the  one  obtained  experimentally. 

3.  Obtain  a  state-space  realization  of  the  low-order  transfer  function  obtained  in  step  2. 

We  assume  that  the  uncontrolled  n-th  flutter  mode  behaves  like  a  lightly  damped  harmonic  oscil¬ 
lator  with  individual  blades  moving  in  the  stationary  frame  according  to  the  formula  (2.44).  Thus, 
we  expect  the  mode  to  have  a  significant  response  to  forcing  only  at  a  narrow  band  of  frequencies  of 
interest  around  the  mode’s  natural  frequency  uns  :=  wnr  +  nwr.  The  control  goal  is  to  add  damping 
to  the  mode  by  a  feedback  control  only  at  this  narrow  band  of  frequencies.  Therefore,  it  is  sufficient 
to  have  an  approximate  low  order  model  describing  the  dynamics  of  the  n-th  mode  at  this  narrow 
frequency  range.  Even  if  the  frequency  response  of  the  n-th  flutter  mode  were  that  of  a  low  pass, 
rather  than  a  band  pass  filter  and  the  actuator  dynamics  cannot  be  neglected  over  a  wide  band  of 
frequencies,  so  that  a  narrow  band  frequency  model  will  not  be  accurate  at  low  frequencies,  the  in¬ 
accuracy  of -the  model  will  not  significantly  impact  control  performance.  The  controllers  will  have  a 
band  pass  characteristic,  so  that  the  unmodelled  dynamics  at  both  low  and  high  frequencies  will  not 
be  destabilized. 

The  transfer  function  between  the  n-th  SFC’s  of  the  forcing  function  and  the  corresponding  blade 
deflection  response  is  defined  by 


an  (j^)  _  An 
Un{juj)  Un 


(2.50) 


Both  An  and  </>n  are,  in  general,  functions  of  the  frequency  u. 

To  obtain  the  transfer  function  Gn(juj)  from  a  sine  sweep  experiment,  one  has  to  access  the 
function  3n(£).  To  obtain  an  approximation  to  5n(£)  one  would  have  to  simultaneously  measure  the 
blade  displacement  an(t,  6)  at  some  finite  number  of  angles  around  the  annulus  and  use  a  discrete 
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approximation  of  the  integral  defining  the  spatial  Fourier  transform.  A  reasonable  approximation 
would  require  at  least  2n  +  1  blade  displacement  sensors  around  the  annulus. 

However,  even  with  one  sensor  one  can  measure  the  transfer  function  Gn (ju)  because  of  the 
following  simple  observation.  Assume  that  we  have  a  blade  displacement  sensor  at  some  angle  0y  at 
the  stationary  frame.  The  measured  output  function  ynov(t)  :==  an{t,9y)  is  given  by  the  equation 
2.49.  Assume  also  that  we  measure  the  value  of  the  actuation  function  un(£,0)  at  a  fixed  angle  0U. 
bet  un0u(t)  :=  UnfaOu).  Note  that  yney(t)  =  Ancos(u)0t  -  nOy  +  <j>n)  and  uneu(t)  =  Uncos(u)0t  -f 
(t>nu  -  n9u)  have  relative  phase  shift  of  <f>n  -  (j)nu  -  n(Qy  -  6U).  Hence,  the  measured  transfer  function 
between  them  during  a  sine  sweep  experiment  is  Gn$ugv{juj)  :=  |n*y =  An .ej(<f>n-<l>nu-n(0y-9u))  _ 

e~in(9v~0u)Gn(ju>).  Therefore,  Gn(ju)  can  be  obtained  from  Gn(juj)  =  e^9y~9u^ Gneuey(jcj) . 

One  can  observe  that,  except  for  the  case  n  =  0,  any  description  of  transfer  functions  Gn(ju) 
as  a  rational  function  of  jcv  valid  in  a  wide  frequency  band  must  have  complex  rather  than  real 
coefficients.  To  see  that,  note  that  a  transfer  function  G(jcj)  with  real  coefficients  has  the  property 
G(—juj)  =  G(ju )*,  i.e.,  it  has  a  Nyquist  diagram  symmeric  with  respect  to  the  real  axis.  We  know 
from  experiments  that  the  response  of  the  n-th  flutter  mode  to  the  forward  or  backward  travelling 
forcing  wave  with  the  same  temporal  frequency  is  not  symmetric.  Thus,  for  n  ±  0,  one  has  Gn{—ju))  ^ 
Gn (jw)*.  However,  we  do  expect  the  reponse  to  be  symmetric  for  n  =  0,  so  that  we  have  Go(juj)  = 
Go{—ju>)*.  Therefore,  we  expect  the  transfer  function  Go(ju>)  considered  as  a  rational  function  of  ju 
to  have  real  coefficients.  Because  of  this  difference  between  the  cases  n  ^  0  and  n  =  0,  we  are  going 
to  derive  the  corresponding  models  separately. 

A  low  order,  narrow  band  model  for  the  transfer  function  Gn(ju>)  between  the  n-th  SFC  of  the 
forcing  function  given  by  (2.47)  and  the  corresponding  output  function  given  by  (2.49)  for  n  ^  0  is  a 
first  order  transfer  function  with  complex  coefficients 


Gn(jv)  = 


bnR  +  jbnl 


Cn+j{v-Uns)' 

A  complex- valued  state-space  realization  of  this  transfer  function  is 


(2.51) 


Qfn(t)  (  Cn  H"  j^ns)&n(j')  H"  (finR  jbnl)un{t) •  (2.52) 

Note  that  both  5 n(t)  and  un(t)  are  complex  valued  functions  of  time.  Observe  also  that  the  unforced 
response  of  (2.52)  is  5 n(t)  —  e(”^n+-?CJna)tan(0),  which  agrees  with  postulated  unforced  evolution  of  the 
n-th  flutter  mode  given  by  (2.44). 

Let  us  emphasize  again  that  the  simple  transfer  function  model  (2.51)  and  its  state-space  realization 
(2.52)  are  valid  only  for  a  narrow  range  of  frequencies  around  the  flutter  frequency  uns.  The  actuator 
characteristic  over  that  frequency  range  is  simply  represented  by  the  magnitude  and  phase  of  the 
n-th  mode  of  the  actuator  disk  at  the  flutter  frequency  and  incorporated  into  the  complex  number 
bnR  +  jbni.  This  approximation  is  reasonable,  as  long  as  the  actuator  frequency  response  does  not 
change  significantly  over  the  frequency  interval  of  interest  and  a  feedback  controller  characteristic  will 
be  that  of  a  sufficiently  narrow  band-pass  filter.  If  this  is  not  the  case,  the  actuator  dynamics  should 
be  incorporated  in  the  model. 

An  equivalent  description  to  (2.51)  is  possible  with  a  real- valued  model  of  real  dimension  two.  In 
the  sequel  the  subscripts  (-)r  and  (•)/  will  denote  the  real  and  imaginary  part  of  a  complex  number. 
One  can  easily  check  [42]  that  the  real  and  imaginary  part  of  the  SFC’s  of  blade  displacement  and 
forcing  function  satisfy  the  following  set  of  two  differential  equations 


anR(t) 

~Cn  1 

&nR(t) 

Unlit) 

^ns  ( 

n  J 

Unl{t) 

bnR  ~bnl 
bnl  bnR 

=?>  f  * 

^  ft 

i _ _ _ _ i 

(2.53) 
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The  corresponding  transfer  function  desription  is 


a„R(jw) 

_  otnI{ju}) 

Gnr(ju)  -Gni(juj) 

Gniijbj)  Gnr(ju)) 


UnR(ju) 

Unl(ju) 


One  can  verify  that 


Gnr{ju))  — 
Gni{jiv)  = 


KrU^  +  Cn)  ~  bul^ns 
(jV  +  Cn)2  +  U)%s 
~^n/(j^  H~  Cn)  ~~  bnRLOns 
[jv  +  Cn)2  +  uls 


(2.54) 

(2.55) 

(2.56) 


and 


Gn(ju)  =  Gnr(juj)  4-  jGni{juj).  (2.57) 

Assume  that  a  blade  displacement  sensor  is  located  at  some  angle  9  at  the  stationary  frame.  The 
measured  output  function  yn$(t)  :==  orn(£,  9)  can  be  expressed  in  terms  of  the  real  and  imaginary  parts 
of  the  SFC  of  an(t,0)  via  the  inverse  spatial  Fourier  transform  (2.45)  as  yne(t)  =  Re(an(t)e~^n$)  = 
Re(anR{t)  +joLni{t)){ cos(n0)  -  j  sin (n6))  —  cos(n0)5n;*(t)  +  sin(n0)5n/(t).  Let 


^n(^)  • — 


5  nR{t) 
&nl(t) 


j vn(t)  :~ 


UnR(t) 

Unl{t) 


(2.58) 


An  := 


,  Bn  : — 


bnR  —bnT 
bnl  bnR 


Cn  ~^ns 
Wns  ”Cn 

Cno  '=  [cos(n0)  sin(n0)]. 

The  state  and  the  output  equation  for  the  n-th  nodal  flutter  mode  can  be  concisely  written  as 


(2.59) 


2*nW  —  Anxn(t)  4*  Bnvn(t) 
yn0(t)  =  CneXn(t). 


(2.60) 


If  there  is  only  one  sensor  at  some  fixed  angle  0 ,  we  will  skip  the  subscript  0  in  the  description  of 
yne(t)  and  CnQ . 

Observe  that  all  the  quantities  in  the  equation  (2.60)  are  real.  One  can  identify  the  parameters 
in  the  model  using  travelling  wave  excitation,  as  described  in  the  previous  section.  Alternatively,  one 
can  exploit  the  skew-symmetric  structure  of  the  matrices  An  and  Bn  and  use  only  one  of  the  inputs  of 
vn(t)  for  excitation.  This  amounts  to  forcing  the  system  with  a  standing  wave,  rather  than  travelling 
wave  pattern. 

Now  we  propose  a  real-valued  model  for  control  of  the  0-th  nodal  flutter  of  dimension  two.  Assume 
that  a  blade  displacement  sensor  is  located  at  some  angle  9  at  the  stationary  frame.  The  measured 
output  function  is  yoe{t)  :=  ao(t,0)  —  50(t).  Let 


.  a0  (juj) 

Go{ju)  :=  , 

uq{3v) 


UQ  {ju)  C/q 


(2.61) 


A  simplest  model  for  Go(jco)  with  real  coeffcients  that  exhibits  a  behavior  of  a  lightly  damped  oscillator 

for  some  real  bi,  bo,  Co>  and  u> o-  The  corresponding  state-space  description  (in  the  observer  canonical 
form)  is 

xo(t)  =  Aoxo(t)  +  Bov0(t) 
yoe(t)  =  C0gx0(t), 


(2.63) 
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where 


(2.64) 


v0(t)  :=  u0(t), 

4>:=  ~^2  l  ,B0:=  ^  ,Coe  :=  [10],  (2.65) 

Let  consider  a  finite  number  kj  of  flutter  modes  with  nodal  numbers  n  i,  n2,  . . Assume  that 
as  the  measurement  outputs  we  use  m  blade  displacement  sensors  located  at  the  angles  02,  •  • 
respectively.  We  assume  that  the  blade  displacement  ye  measured  at  some  angle  0  at  the  stationary 
frame  is  the  sum  of  the  diplacements  due  to  particular  flutter  modes: 


yo(t)  =  Efc=l  Olnk(t,  9). 


(2.66) 


We  can  write  down  the 
flutter  modes: 


following  state-space  model  describing  the 

x(t)  =  Ax(t)  -f  Bv(t) 
y{t)  =  Cx(t), 


dynamics  of  the  kj  most  active 


(2.67) 


where  x{t)  ~  [xni  (t) . . .  xnkf  (t)]r, 

v(t)  :=  [uni  (<)...  unkf  (i)]T,  y(t)  [y8l  (t) . . .  y6m  (t)]T,  A  and  B  are  block  diagonal  matrices  containing 
Anj  and  Bnj  blocks,  respectively,  and  C  is  a  matrix  composed  of  Cne  blocks. 

The  dimension  of  the  output  variable  y(t)  is  equal  to  m,  which  is  the  number  of  blade  displacement 
sensors  (e.g.,  eddy  current  sensors)  used  for  measurement. 

One  may  be  tempted  to  place  many  sensors  to  make  the  C  matrix  invertible  and  use  a  full- 
state  static  feedback  to  arbitrarily  place  the  damping  of  the  flutter  modes.  This  strategy  might  be 
succesfull  for  flutter  modes  with  the  nodal  number  ^  0  if  the  variations  of  the  actuator  dynamics 
with  frequency  can  be  neglected.  However,  note  that  C  is  never  invertible  if  one  includes  the  0-th 
nodal  flutter  dynamics,  as  Coe  =  [1  0]  for  all  0,  and  hence  C  has  a  column  of  zeros.  Moreover,  a 
strong  output  noise  component,  which  includes  all  unmodelled  sources  of  blade  displacement,  such  as 
periodic  forcing  due  to  rotor  and  blades  assymetry,  neglected  flutter  modes,  rotating  stall  dynamics, 
an  inlet  distorsion,  etc.,  would  make  reconstructing  the  state  of  the  flutter  modes  by  inverting  the  C 
matrix  problematic. 

To  circumvent  the  problems  with  direct  inversion  of  the  C  matrix  and  the  output  noise,  and  at 
the  same  time  reduce  the  number  of  blade  displacement  sensors,  we  are  going  to  reconstruct  the  state 
of  the  system  using  an  observer.  As  we  will  see,  in  principle,  just  one  blade  displacement  sensor  is 
sufficient  for  this  purpose. 


2.5.4  Model  of  the  disturbances 

We  are  going  to  augment  the  state-space  model  (2.67)  by  adding  some  noise  sources  in  the  state  and 
output  equations.  We  assume  that  the  system  is  described  by  the  equation 


x(t)  =  Ax(t)  Bv(t)  -f  ex(t) 
y(t)  =  Cx(t)+ey{t), 


(2.68) 


where  ex(t)  is  an  unknown  disturbance  driving  the  state  of  the  system,  and  ey(t)  is  an  unknown 
output  disturbance.  We  assume  that  the  state  disturbance  ex(t)  has  a  strong  periodic  component 
with  the  rotor  frequency  due  to  rotor  and  blades  assymetry  and  some  random  fluctuations  due  to  inlet 
distorsion  that  excite  the  kf  flutter  modes  of  interest.  The  output  disturbance  ey(t)  might  include  all 
unmodelled  sources  of  blade  displacement  such  as  neglected  flutter  modes,  rotating  stall  dynamics, 
a  strong  constant  component  (so-called  DC-component)  due  to  sensor  bias  and  rotor  assymetry,  a 
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periodic  component  due  to  differences  between  particular  blades,  60 Hz  electric  noise  and  its  harmonics 
from  electric  components  and  possibly  some  other  measurement  noise. 

We  expect  the  constant  and  periodic  components  of  the  disturbances  to  be  dominant.  We  are 
going  to  extend  our  model  to  incorporate  these  components  and  reconstruct  them  using  an  observer. 
We  assume  that  the  actual  blade  displacement  yae{t)  measured  by  a  sensor  located  at  an  angle  9  on 
the  stator  is  given  by  yae{t)  =  y$(t)  4-  ydc{&)  +  yP{vrt  -  0)  +  ywn{t,9),  where  ye(t)  is  the  sum  of  the 
blade  displacement  components  due  to  the  dynamics  of  the  selected  most  active  kf  flutter  modes  of 
interests  given  by  (2.66),  ydc{9)  is  an  unknown  constant  bias  component,  yp{wrt  —  6)  represents  a 
steady  unknown  periodic  motion,  and  yun{t ,  6)  is  the  broad  band  component  of  the  measured  output. 

Note  that  the  forcing  of  the  blades  due  to  rotor  and  casing  assymetry  will  cause  a  different  constant 
blade  diplacement  at  a  different  angular  location  at  the  stationary  frame,  even  the  the  magnitude  of 
the  angular  motion  of  the  blades  is  the  same.  Also,  a  constant  bias  signal  may  be  different  fo  each 
sensor.  Therefore,  we  do  not  expect  to  have  the  same  constant  component  a<fc(0)  at  each  output. 

The  unknown  constant  component  of  the  fc-th  output  yok(t)  will  be  modelled  as  the  state  x^dcit)  of 
an  integrator  with  an  unknown  initial  condition  added  to  the  system  (2.68)  xkdc(t)  =  AdcXdc(t),  Vdk  (t)  = 
CdcXkdc(t),  where  :=  [0],  Qc  :=  [1]. 

Particular  blades  on  the  rotor  are  slightly  different.  Moreover,  the  spacing  between  the  blades  may 
vary  around  the  rotor.  Thus,  if  the  diplacement  is  measured  by  a  difference  of  expected  (assuming 
no  blade  motion  due  to  flutter)  and  actual  time  of  arrival  of  the  blade  (using  eddy  current  or  optical 
senors)  and  the  time  of  expected  arrival  is  calculated  assuming  a  uniform  spacing  between  the  blades, 
a  measurement  of  an  angular  blade  displacement  will  have  an  error.  In  both  cases  mentioned  above 
this  error  will  have  a  form  of  a  wave  travelling  forward  with  the  speed  equal  to  the  rotor  frequency 
ay.  This  justifies  the  choice  of  the  form  of  the  periodic  disturbance  yp(urt  —  0). 

A  constant  component  of  yp{u)rt  —  6)  can  be  incorporated  into  the  constant  disturbance  c*dc(0). 
Therefore,  without  loss  of  generality  one  can  assume  that  for  each  fixed  t  the  average  value  of  yp(wrt— 9) 
over  9  is  zero.  For  each  n  >  0  one  can  compute  the  spatial  Fourier  coefficient  ypn(t)  of  yp(oyt  —  9)  and 
then  represent  the  function  yp(u)rt  —  9)  using  the  inverse  spatial  Fourier  transform  (in  other  words, 
represent  yp(u)rt  —  9)  as  Fourier  series).  After  truncating  the  series  at  some  n  =  kp  one  obtains  the 
formula  yp(wrt  -  9)  =  £n=i  Re(ypn(t)e~jne)  =  £n=l  cos(n0)(a„  cos(na;rt)  4-  bn  sin(nayt)) 

+  sin(n0)(an  sin(nuyt)  —  bn  cos  (nay  t)),  for  some  constants  on,  6„,  n  =  1, . . . ,  kp.  Observe  that  for  fixed 
0,  the  function  yp(oy£  ~  0)  can  be  treated  as  an  output  of  a  system  of  kp  uncoupled  oscillators  with 
frequencies  equal  to  multiples  of  the  rotor  frequency  ay.  Therefore,  the  periodic  disturbance  yp{wrt— 9) 
will  be  modelled  by  adding  states  xie(t),  X2e{t), . . .  ,xicpe{t)  of  kp  oscillators  to  the  system  (2.68).  Let 
tjfce  :=  kujr  and 

Ake  :=  °  J ke  ,Cn*  =  [cos(n0)  sin(n0)].  (2.69) 

U)ke  o 

We  assume  that  the  state  of  the  oscillator  xke(t)  satisfies  the  equation 

^fce(^)  =  Akexke{£))  (2.70) 

with  an  unknown  initial  condition.  The  function  yp(ujrt  —  9)  can  be  expressed  as 

kp 

yp(L>Tt  -0)=  Y,  Cnf>Xne(t)  (2.71) 

n—l 

We  assume  that  yWn{ti9),  the  unmodeled  part  of  the  output  function,  is  a  low  level  broad  band 
noise. 

We  add  the  disturbance  states  x\ rfc(t),  •  •  • ,  awc(t)  and  xie(t), . . . ,  Xkpe{t)  to  the  state  x(t)  and  form 
an  augmented  state  variable  xa(t).  We  can  write  down  the  following  augmented  state-space  model 
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describing  the  dynamics  of  the  kf  most  active  flutter  modes  and  the  effect  of  disturbances  on  the 
output  and  state: 

xa(t)  =  Aaxa(t)  +  Bav(t)  +  eaXu{t)  ( 9  79v 

Va{t)  =  Caxa(t)  +  eayu(t),  1  ' 

where  eaxw(t)  and  eayu(t)  denote  the  broad-band  components  of  the  state  and  output  disturbances 
(we  skip  the  details  of  straightforward  contruction  of  the  augmented  system  matrices). 

Note  that  the  system  (2.72)  is  not  controllable,  as  we  do  not  assume  that  we  have  any  actuation 
authority  over  the  part  of  the  system  that  models  the  disturbances.  Therefore,  we  can  only  attempt 
to  place  the  eigenvalues  of  the  system  corresponding  to  the  flutter  modes.  If  the  full  state  of  (2.72) 
were  available  for  measurement,  we  would  use  a  feedback  law  of  the  form  v  :=  —Kaxay  where  Ka  is  of 
the  form  Ka  —  [K  0]  and  K  is  chosen  so  that  the  eigenvalues  of  the  matrix  A  —  BK  are  placed  in  the 
desired  position.  This  can  be  achieved,  as  it  follows  from  the  Hautus  criterion  that  the  pair  ( A ,  B)  is 
controllable  since  the  flutter  modes  are  separated  in  frequency.  Note  that  the  matrices  A  and  B  are 
block-diagonal.  Hence,  one  can  use  a  block  diagonal  feedback  gain  matrix  K  and  separately  place  the 
eigenvalues  of  particular  flutter  modes.  As  we  noticed,  the  state  of  the  flutter  modes  is  not  directly 
available.  However,  one  can  use  an  observer  state  instead  of  the  actual  state  for  feedback. 


2.5.5  Observer-based  control  of  flutter 

The  main  reason  for  incorporating  the  unknown  disturbances  into  the  state  of  the  system  is  to  re¬ 
construct  them  along  with  the  state  of  the  flutter  dynamics  and  filter  out  in  this  way  the  state  of 
the  flutter  modes  of  interest.  This  state  can  be  used  for  designing  a  full  state  feedback  and  achieve  a 
prescribed  level  of  damping  of  the  flutter  modes. 

To  verify  if  the  augmented  system  (2.72)  is  observable  one  can  use  the  Hautus  criterion  of  observ¬ 
ability  [29],  The  latter  says  that  (2.72)  if  and  only  if  the  matrix 


A  I-Aa 
Ca 


(2.73) 


has  a  full  column  rank  VA  E  a (Aa),  where  cr(-)  stands  for  the  spectrum  of  a  matrix.  Note  that  the 
matrix  Aa  is  block-diagonal.  Thus,  its  spectrum  is  the  collection  of  spectra  of  the  matrices  on  the 
diagonal.  One  can  easily  see  that  as  long  as  the  flutter  frequencies  do  not  coincide  with  multiples  of 
the  rotor  frequency,  the  system  is  observable  even  from  a  single  sensor. 

In  the  sequel  we  assume  that  the  multiples  of  the  rotor  frequency  are  separated  from  the  frequencies 
of  the  flutter  modes,  so  that  the  pair  (Cay  Aa)  is  observable. 

An  observer  used  to  reconstruct  flutter  modes  [29]  has  the  form 

xa(t)  =  Aaxa{t)  +Bava(t)  +  L{ya{t)  -ya{t)) 

y{t)  =  CaXait),  ^  ' 

with  the  observer  gain  matrix  L  chosen  so  that  the  observer  error  converges  to  zero.  Let  x(t)  := 
xa{t)  —  xa(t)  denote  the  error  of  observation.  The  error  dynamics  is  given  by 

Xa (^)  =  (Aa  LCa):ra(t)  LeayU(t)  -h  eamt(£)*  (2.75) 

To  guarantee  the  convergence  of  the  observer  error  to  zero,  the  matrix  L  has  been  chosen  such  that  the 
matrix  Aa  —  LCa  has  all  eigenvalues  with  negative  real  part.  Since  the  pair  ( Ca ,  Aa)  is  observable,  the 
eigenvalues  of  Aa  —  LCa  can  be  assigned  arbitrarily.  By  choosing  the  eigenvalues  far  enough  from  the 
imaginary  axis  one  can  make  the  observer  error  decaying  arbitrarily  fast.  However,  this  requires  high 
values  of  the  entries  of  the  observer  gain  matrix  L.  Since  the  unmodelled  part  of  the  measurement 
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noise  eayu(t )  affects  the  error  dynamics  (2.75)  through  the  matrix  L,  a  high  gain  observer  is  going  to 
have  a  significant  response  to  the  measurement  noise.  Therefore,  one  has  to  choose  the  eigenvalues  of 
the  matrix  Aa  —  LCa  carefully,  balancing  the  rate  of  the  observer  error  decay  with  the  sensitivity  to 
output  noise. 

One  can  design  a  feedback  controller  for  the  system  (2.72)  in  a  form  of  a  full-state  feedback 
va(t)  =  —Kaxa(t),  where  Ka  =  [K  0]  is  chosen  so  that  the  eigenvalues  of  the  matrix  A  —  BK 
corresponding  to  the  flutter  modes  have  a  prescribed  level  of  damping  [29].  Since  the  full  state  of 
(2.72)  is  not  accessible  for  a  direct  measurement,  one  uses  the  state  of  the  observer  (2.74)  in  place 
of  the  state  of  (2.72)  in  the  feedback  law.  i.e.,  va(t)  =  —Kaxa(t).  The  state  equations  of  the  system 
consisting  of  an  interconnection  of  the  system  (2.72)  and  the  observer  (2.74)  are 


xa{t) 

-BaKa 

Xa{t) 

xa(t) 

“  1  J 

S' 
_ 1 

-  BaKa  -  LCa 

_  Xa(t) 

/ 

0 


eaxu(t)  + 


+ 


(2.76) 


The  Separation  Principle  says  that  the  set  of  eigenvalues  of  the  overall  system  (2.76)  is  a  union  of  the 
eigenvalues  of  the  observer  error  matrix  Aa  —  LCa  and  the  eigenvalues  of  the  matrix  Aa  —  BaKa.  In 
turn,  the  set  of  eigenvalues  of  Aa  —  BaKa  is  the  union  of  the  eigenvalues  of  A  —  BK  and  the  eigenvalues 
of  the  modelled  disturbance  modes. 


2.5.6  Flutter  control  experiments 

Experiments  with  observer-based  flutter  control  on  17" -fan  test  rig  shown  in  Figure  2.13  .  The  rig  has 


Figure  2.13:  17”  fan  experimental  rig 

one  fan  stage  with  sixteen  blades  on  a  rotor  powered  with  an  external  air  turbine  and  one  row  of  exit 
guide  vanes  on  a  stator.  Variable  fan  exit  area  is  controlled  via  translating  throttle  plug.  By  reducing 
the  throttle  area  one  increases  the  incidence  angle  of  the  air  approaching  the  blades  and  thus  reduces 
damping  of  stall  flutter  modes.  The  control  system  is  shown  in  Figure  2.14.  We  used  two  or  one  eddy 
current  sensors  to  measure  the  blade  displacement  in  the  stationary  frame  and  ten  speakers  mounted 
in  cavities  located  around  the  fan  casing  as  actuators.  All  control  experiments  reported  here  done  at 
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Off  Blade  Active  Flutter  Control  Schematic 

.  Digital  Counting  Circuit 


Figure  2.14:  Flutter  control  system  schematics 


the  rotor  frequency  150  Hz  (9000  RPM).  Most  expreriments  were  done  at  the  most  stable  operating 
point  at  wide  open  throttle  and  some  closer  to  flutter  boundary. 

The  goal  was  to  demonstrate  that  we  can  add  significant  amount  of  damping  to  the  lightly  damped 
flutter  modes.  The  goal  was  achieved  for  the  least  stable  fan  flutter  modes  with  nodal  diameters  0,  1, 
and  2. 

We  identified  the  parameters  of  the  real  valued  2-nd  order  state-space  models  for  the  flutter  modes 
0,  1,  and  2  from  a  sine  sweep  experiment.  For  nodal  diameter  0  we  assumed  model  given  by  equation 
(2.63).  For  the  nodal  diameters  other  than  0  we  assumed  model  given  by  equation  (2.60).  We  used 
a  Matlab  script  that  first  fitted  a  2-nd  order  transfer  function  (with  two  poles  and  one  zero)  to  the 
experimentally  obtained  transfer  fucntion  from  a  spectral  analyzer  and  then  obtained  an  appropriate 
state-space  realization.  For  nodal  diameters  n  =  1  and  n  =  2  we  used  forcing  with  a  standing  wave 
cosine  pattern,  i.e.,  the  input  unn(t).  In  this  way  we  identified  the  matrices  Cuq  (for  two  values  of 
0),  Anj  and  the  first  column  of  the  matrix  £n.  The  matrix  An  was  forced  to  be  in  the  (modal)  form 
(2.59).  The  second  column  of  Bn  was  obtained  from  the  first  one. 

Figures  2.15,  2.16  and  2.17  show  the  Bode  plots  from  a  sine-sweep  experiment  and  the  correspond¬ 
ing  2-nd  order  fits  for  flutter  modes  with  nodal  diameters  n  =  0,  n  =  1,  and  n  =  2.  The  fits  and  models 
are  obtained  independently  for  each  of  the  eddy  current  sensors  (denoted  DTI  and  DT2,  respectively). 

Note  that  the  2-nd  order  transfer  function  fits  to  experimental  data  are  good  near  the  flutter 
frequency  (i.e.,  the  frequency  at  the  resonant  peak  of  the  magnitude),  but  there  is  an  increasing  phase 
error  away  from  the  flutter  frequency.  Ass  we  will  see,  this  phase  mismatch  represents  unmodeled 
non-minumum  phase  dynamics  that  will  somewhat  limit  performance  of  the  control  scheme  based  on 
the  2-nd  order  model. 

Several  observer-based  schemes  we  tried  in  flutter  control.  Some  of  them  used  both  sensors  and 
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Figure  2.15:  Bode  plots  for  flutter  mode  n  =  0.  Experimental  data  and  2-nd  order  fits. 


some  of  them  used  just  one  sensor.  Number  of  flutter  and  disturbance  modes  also  varied.  All  observers 
had  a  constant  disturbance  states  (denoted  dc)  and  from  one  to  four  pairs  of  the  states  modelling 
periodic  disturbances  at  the  multiples  of  the  rotor  frequency  (denoted  le,  2e,  etc.).  The  available 
processor  speed  limited  the  number  of  observer  states  that  could  be  used  for  control. 

Each  observer  was  designed  in  the  form  (2.74),  with  the  gain  matrix  L  chosen  by  specifying  the 
damping  of  the  observer  poles.  The  specified  observer  damping  was  about  twice  the  value  of  the  desired 
damping  of  the  flutter  mode  of  interest,  except  for  the  dc  states,  which  had  much  higher  damping. 
The  control  law  was  a  linear  feedback  from  the  observer  state.  Only  the  (estimated)  flutter  state  was 
used  for  control.  The  feedback  was  designed  to  augment  the  damping  of  the  flutter  mode  of  interest 
to  the  desired  level  without  changing  its  frequency. 

The  amount  of  damping  added  to  particular  flutter  modes  was  obtained  from  a  sine-sweep  ex¬ 
periment.  Bode  plots  for  open  and  closed-loop  systems  were  obtained.  Since  the  2-nd  order  transfer 
function  fits  were  not  acurate  enough  (see  Figures  2.15,  2.16  and  2.17  )  to  get  a  good  estimate  of 
damping,  4-th  order  transfer  functions  were  fitted  instead.  The  4-th  order  fits  were  very  good  for  0-th 
and  1-st  nodal  diameter  flutter  modes,  and  reasonable  for  the  2-nd  nodal  diameter  flutter  mode. 

The  0-th  flutter  mode  was  the  least  stable  mode  on  17-inch  fan  rig.  At  the  most  stable  rig  operating 
point  the  open-loop  logarithmic  decrement  was  Jo  =  0.018.  The  control  scheme  that  achieved  the 
biggest  damping  increase  used  only  one  DT  sensor  and  was  based  on  model  with  7  states:  2-nd  order 
models  of  the  0-th  and  1-st  flutter  modes,  1-st  order  model  of  the  dc  bias  of  DTI  sensor,  and  a  2-nd 
order  model  of  2e  signal.  The  closed-loop  logarithmic  decrement  achieved  was  Jo  =0.21. 

Figure  2.18  shows  open  and  closed-loop  Bode  plots  for  the  plant  and  the  4-th  order  transfer  function 
fits. 

Figure  2.19  shows  the  spectral  content  of  the  time  traces  of  one  actuator,  strain  gauge  on  a  blade, 
and  DTI  (time  difference  between  actual  and  expected  blade  arrival)  signal.  The  frequency  of  the 
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Figure  2.16:  Bode  plots  for  flutter  mode  n  =  L  Experimental  data  and  2-nd  order  fits. 


0-th  flutter  mode  (about  273  Hz)  is  shown  by  a  circle.  The  first  three  multiples  of  the  rotor  frequency 
(150  Hz)  are  marked  by  crosses.  Note  that  the  actuator  has  a  strong  component  at  the  (unmodeled) 
le  frequency,  but  no  component  at  the  (modeled)  2e  frequency,  even  though  there  is  a  strong  2e  signal 
in  the  DTI  signal.  This  illustrates  the  benefits  of  estimating  the  periodic  disturbances  at  multiples 
of  rotor  frequency  using  an  observer.  Essentially,  the  periodic  disturbance  observer  provides  a  notch 
filter  in  the  transfer  functions  from  blade  arrival  time  difference  to  the  speaker  command.  Adding 
more  models  of  the  periodic  disturbances  with  muliples  of  the  rotor  frequency  to  the  observer  reduces 
the  actuator  energy  waisted  for  response  in  that  frequencies. 

The  1-st  flutter  mode  was  the  second  least  stable  mode  on  17-inch  fan  rig.  At  most  stable  rig 
operating  point  the  open-loop  logarithmic  decrement  was  Si  =  0.028.  The  control  scheme  that  achieved 
the  biggest  damping  increase  used  both  DTI  and  DT2  sensors  and  was  based  on  model  that  had  8 
states:  2-nd  order  models  of  0-th  and  1-st  flutter  modes,  dc  estimates  for  DTI  and  DT2  sensors,  and 
a  2-nd  order  model  of  2e  signal.  The  closed-loop  logarithmic  decrement  was  8\  =  .153.  Figure  2.20 
shows  open  and  closed-loop  Bode  plots  for  the  plant  and  the  4-th  order  transfer  function  fits. 

The  2-nd  flutter  mode  was  the  third  least  stable  mode  on  17-inch  fan  rig.  At  the  most  stable  rig 
operating  point  the  open-loop  logarithmic  decrement  was  <5i  =  0.037.  We  used  both  DTI  and  DT2 
sensors.  The  model  used  for  observer  design  had  10  states:  2-nd  order  model  of  2-nd  flutter  mode, 
dc  estimates  for  DTI  and  DT2  sensors,  and  2-nd  order  models  of  3e,  4e,  and  5e  disturbance  signals. 
The  closed-loop  logarithmic  decrement  was  =  .174.  Figure  2.21  shows  open  and  closed-loop  Bode 
plots  for  the  plant  and  the  4-th  order  transfer  function  fits. 

Figure  2.22  summarizes  achieved  damping  augmentation. 

We  tested  an  observer-based  controllers  designed  at  the  most  stable  operating  point  at  the  rotor 
frequency  of  9000  RPM  for  the  throttle  position  very  close  to  the  stability  boundary  of  the  n  =  0 
flutter  mode.  Figure  11  shows  the  open  and  closed-loop  Bode  plots  and  damping  estimates.  The 
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Figure  2.17:  Bode  plots  for  flutter  mode  n  =  2.  Experimental  data  and  2-nd  order  fits. 


controller  worked  well.  This  was  expected,  as  the  phase  of  the  frequency  response  of  the  0-th  flutter 
mode  did  not  change  much  as  one  closed  the  throttle. 

At  the  same  operating  condition  (close  to  the  0-th  flutter  mode  stability  boundary)  we  conducted 
a  control  off/on/off  experiment.  Figure  A2.24  shows  10-second  time  traces  of  one  of  the  actuators,  a 
strain  gauge  on  a  blade,  and  DTI  signal.  One  can  see  a  significant  decrease  in  the  blade  stress  and 
blade  tip  movement  due  to  0-th  flutter  mode  when  the  feedback  loop  is  closed.  Figure  2.25  shows  the 
spectral  content  of  these  time  traces  when  control  was  off  and  on,  respectively.  Note  the  disappearance 
of  a  sharp  peak  at  273  Hz  at  both  the  strain  gauge  and  the  blade  arrival  time  difference  signal  spectra 
when  control  is  on. 

2.5.7  Experimental  results:  comments  on  model  mismatch  effect 

The  amount  of  damping  augmentation  achieved  for  the  flutter  modes  with  nodal  diameters  0,  1,  and 
2,  was  significant.  A  prediction  based  on  the  2-nd  order  models  was  that  even  more  damping  could 
have  been  achieved  simply  by  increasing  the  controller  gain.  However,  in  experiment  we  noticed  that 
there  was  an  upper  bound  on  the  damping  augmentation.  As  we  were  increasing  the  gain,  the  plant 
frequency  response  for  all  flutter  modes  consistently  showed  a  magnitude  plot  with  two  resonant  peaks, 
characteristic  to  a  4-th  order  system.  Indeed,  the  4-th  order  transfer  function  fits  (with  4  poles  and 
3  zeros)  to  the  experimental  transfer  functions  were  excellent  over  40  —  50  Hz  frequency  band  around 
the  resonant  peak  for  both  open-loop  and  closed-loop  plants  for  all  values  of  the  controller  gains  that 
we  tried.  For  each  open-loop  flutter  mode  the  four-pole  pattern  was  similar.  There  were  two  complex 
conjugate  pairs  of  poles  with  a  similar  frequency.  One  pair  was  located  near  the  imaginary  axis,  more 
or  less  at  the  location  of  the  poles  of  the  second  order  fit.  The  other  pair  of  poles  was  located  further 
away.  Our  interpretation  was  that  the  less  damped  pair  of  poles  represented  the  flutter  dynamics, 
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Figure  2.18:  Open-loop  and  closed-loop  Bode  plots  for  flutter  mode  n  =  0.  Experimental  data  and 
4-th  order  fits.  Sq :  open-loop  0.018,  closed-loop  0.21. 


while  the  more  damped  extra  pair  of  poles,  together  with  an  acompanying  pair  of  zeros,  accounted 
for  the  actuator  dynamics  and  all  sources  of  delay  in  the  loop.  It  is  important  to  note  that  we  do 
not  associate  any  direct  physical  meaning  with  the  ’’actuator”  poles  and  zeros.  Rather  than  that,  we 
interpret  them  as  a  narrow  band  bulk  model  of  the  phase  mismatch  between  the  experimental  transfer 
functions  and  the  2-nd  order  fits,  with  the  actuator  phase  rolloff  being  the  major  factor  in  that  roloff. 

As  we  increased  the  controller  gain  the  poles  moved  towards  each  other,  the  flutter  mode  gaining 
stability  and  the  actuator  pole  loosing  stability.  As  a  measure  of  the  damping  of  the  system  we  used 
the  damping  of  the  least  damped  pole  in  the  4-th  order  transfer  function  fit  to  the  experimental 
transfer  function.  The  optimal  damping  was  achieved  as  the  two  pairs  of  poles  nearly  met.  As  we 
kept  increasing  the  gain  from  this  optimal  value,  the  poles  started  to  move  away  from  one  another, 
and  one  of  the  modes  kept  loosing  stability. 

Figure  ??  show  the  experimetally  obtained  open-loop  and  closed-loop  Bode  plots  for  1-st  flut¬ 
ter  mode  and  the  corresponding  4-th  order  fits  for  6  values  of  the  controller  gain  k  (with  k  =  0 
corresponding  to  open- loop  plant).  The  implemeted  modal  control  was  of  the  form 


Vi(t)  := 


«1  R.{t) 
uU(t) 


=  kKi 


xm{t) 

xu{t) 


(2.77) 


where  x \  a  and  xu  were  the  observer  states  corresponding  to  the  1-st  flutter  mode,  K\  was  a  fixed 
feedback  matrix.  The  optimal  value  for  the  gain  k  that  acheieved  largest  damping  augmentation  was 
around  k  =  0.003.  For  the  value  of  gain  higher  than  the  optimal  gain  the  familiar  peak-splitting 
occured.  This  points  out  to  non-minimum  phase  effects  (actuation  and  flow  delays). 

It  was  important  to  verify  that  the  4-th  order  model  would  explain  the  root-locus  behavior  seen 
in  the  experiment.  We  simulated  the  controller  based  on  2-nd  order  model  applied  to  the  4-th  order 
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Figure  2.19:  Control  of  0-th  nodal  diameter  flutter.  FFT  of  time  traces  of  one  actuator,  strain  gauge, 
and  DTI  signal 
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Figure  2.20:  Open-loop  and  closed-loop  Bode  plots  for  flutter  mode  n  —  1.  Experimental  data  and 
4-th  order  fits.  S\:  open-loop  0.028,  closed-loop  0.153. 


model  of  the  1-st  flutter  mode  (the  one  that  showed  the  biggest  phase  mismatch  between  the  2-nd  and 
4-th  order  transfer  function  fits  to  the  experimental  one).  We  plotted  the  corresponding  root  locus  as 
well  as  the  6  flutter-actuator  pairs  of  eigenvalues  obtained  from  the  4-th  order  transfer  function  fits. 
One  can  see  a  qualitative  agreement  between  the  two  root-loci. 
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Figure  2.21:  Open-loop  and  closed-loop  Bode  plots  for  flutter,  n  =  2,  two  eddy  current  sensors 
Experimental  data  and  4-th  order  fits.  5:  open-loop  0.037,  closed-loop  0.174. 
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Figure  2.22:  Flutter  damping  augmentation  achieved  in  experiments 
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Figure  2.24:  Control  off/on/off  experiment,  0-th  nodal  diameter  flutter.  Time  traces  of  one  actuator, 
strain  gauge,  and  DTI  signal. 
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Figure  2.25:  Control  off/on/ofF  experiment,  0-th  nodal  diameter  flutter  close  to  flutter  boundary.  FFT 
of  time  traces  of  one  actuator,  strain  gauge,  and  DTI  signal  for  control  off  (left)  and  control  on  (right). 
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Figure  2.26:  Open-loop  and  closed-loop  Bode  plots  for  1-st  flutter  mode.  Experimental  data  and  4-th 
order  fits  for  6  values  of  the  controller  gain.  Peak-splitting  appears  for  high  gains. 
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Figure  2.27:  Experimental  and  theoretical  root-locus  for  4-th  order  plant  model  and  controller  based 
on  2-order  order  model.  Experimental  closed-loop  poles  for  6  values  of  the  gain  shown  by  circles  and 
crosses.  Theoretical  closed-loop  poles  for  shown  by  crosses. 
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Dear  Dr.  Heise: 

vdr  'f$'* 

Enclosed  is  a  final  report  for  AFOSR  contract  F49620-01-C-0021  "Control  of  Mixing  in  Shear 
Flows".  Please  let  me  know  if  the  report  is  acceptable  in  its  current  form.  Upon  receipt  of  the 
your  acceptance  letter  copies  of  the  report  will  be  sent  to  ACO  and  DTIC. 

Sincerely  yours,  *** 


Andrzej  Banaszuk 

Principal  Research  Engineer 

United  Technologies  Research  Center 


